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Abstract 

Komolgorov,  Arnold  and  Moser  (KAM)  theory  provides  that  orbits  of  satellites  whose 
dynamics  are  representable  by  an  integrable  Hamiltonian  plus  a  small,  real  perturbation 
he  on  tori  in  phase  space  and  remain  upon  the  KAM  tori  for  all  time,  unless  acted  on  by  a 
non-conservative  force.  A  refined  technique  for  constructing  KAM  tori  for  Earth-orbiting 
satellites  is  developed  and  implemented  using  numerically  integrated  orbital  data  for  hy¬ 
pothetical  satellites  and  involving  methods  of  Eourier  analysis  and  spectral  decomposition. 
Definition  of  satellite  formations  on  the  KAM  tori  is  performed  and  analyses  conducted  to 
investigate  both  constellations  with  large  separations  and  clusters  with  small  separations. 
Cluster  formations  with  physical  secular  drift  rates  on  the  order  of  nanometers  to  microm¬ 
eters  per  second  are  obtained.  A  brief  discussion  of  effects  of  non-conservative  forces  (such 
as  atmospheric  drag)  on  KAM  tori  is  given. 
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FORMATION  FLIGHT  OF  EARTH  SATELLITES  ON  KAM  TORI 


I.  Introduction 

There  is  no  question  that  the  advent  and  utilization  of  Earth-orbiting  satellites  in  the 
last  half  century  has  afforded  the  human  race  unprecedented  benefits  and  enhancements  in 
almost  every  aspect  of  its  increasingly  technologically-based  existence.  Among  the  facets 
of  modern  humanity  most  profoundly  affected  by  satellite-based  technology  are  communi¬ 
cations,  navigation,  surveillance  and  reconnaissance,  geographic  analysis  and  deep  space 
exploration. 

1 . 1  Motivation 

In  many  applications  of  satellite  technology,  the  relative  motion  between  two  or  more 
separate  entities  orbiting  about  a  primary  body  is  increasingly  vital.  This  has  become  all 
the  more  apparent  with  the  relatively  recent  focus  on  development  and  implementation  of 
small  satellite  or  microsatellite  formations  -  groups  or  systems  of  smaller,  less  expensive 
satellites  designed  to  perform  or  improve  the  function  of  a  previous  larger,  more  expensive 
satellite  -  for  scientific,  commercial  and  military  purposes. 

Current  methods  of  determining  satellite  relative  motion  and  designing  satellite  con¬ 
stellation  orbits  typically  employ  only  estimates  of  the  true  orbital  dynamics,  to  include 
the  two-body  motion  along  with  the  harmonic  terms  of  Earth’s  gravitational  potential  to 
perhaps  degree/order  m,n  =  4  (see  sections  2.1  and  2.2.1).  While  the  truncation  of  the 
potential  is  sometimes  necessary  due  to  computational  complexity,  especially  during  long 
numerical  integrations,  it  causes  an  increase  in  residual  error  over  time  due  to  the  ignored 
higher-order  gravitational  harmonic  contributions,  except  in  very  specific  cases.  A  definite 
accuracy  advantage  could  obviously  be  realized  in  orbit  design  if  there  exists  a  way  of 
describing  more  fully  and  accurately  the  flight  dynamics  of  the  satellite  about  the  primary. 
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1 . 2  Approach 

In  the  course  of  this  work,  the  author  attempts  to  demonstrate  that,  using  the  con¬ 
cepts  and  results  of  KAM  theory,  it  is  possible  to  significantly  enhance  the  effectiveness  of 
satellite  formation  orbit  design.  This  stems  in  large  part  from  the  fact  that  one  can  achieve 
a  better  representation  of  the  true  dynamics  of  a  satellite’s  orbit  by  using  KAM  theory. 
According  to  the  Kolmogorov,  Arnold  and  Moser  (KAM)  theorem,  the  solution  of  a  near- 
integrable  Hamiltonian  system  will  lie  on  an  invariant  torus  in  the  phase  space  ([1],[8],[15]). 
In  the  case  of  Earth-orbiting  satellites,  if  one  assumes  that  the  only  perturbations  away 
from  the  integrable  two-body  problem  (2BP)  experienced  in  the  orbit  stem  from  Earth’s 
gravitational  field  harmonics  (i.e.  that  nonconservative  perturbations  and  non-primary 
potential  effects  such  as  third-body  interactions  are  zero),  it  has  been  demonstrated  that 
the  satellite’s  motion  may  indeed  lie  on  a  so-called  KAM  torus  for  both  observed  satellite 
data  and  integrated  orbital  data  ([24],  [25], [11]). 

This  work  will  proceed  by  first  demonstrating  the  determination  of  KAM  tori  repre¬ 
senting  various  Earth  satellite  orbits  in  the  presence  of  the  Earth’s  extended  gravitational 
potential.  This  will  be  accomplished  using  techniques  based  upon  those  outlined  in  previ¬ 
ous  works  by  Laskar  in  [9], [10]  and  Wiesel  in  [24], [25]  .  The  application  of  KAM  theory  to 
satellite  formation  flight  will  then  be  investigated  by  exploring  satellite  relative  dynamics 
on  calculated  KAM  tori. 

1.3  Problem  Statement 

Given  the  previously  demonstrated  strong  probability  of  the  existence  of  a  KAM 
torus  upon  which  each  satellite’s  orbit  is  constrained,  the  application  of  KAM  theory  to 
satellite  formation  flight  naturally  arises.  If  the  dynamics  of  earth  orbits  are  indeed  accu¬ 
rately  represented  by  KAM  tori,  the  investigation  of  the  relative  motion  of  two  individual 
satellites  constrained  to  separate  but  proximate  locations  on  the  same  KAM  torus  should 
demonstrate  acceptably  small  secular  drift  between  the  bodies  during  orbital  propagation. 
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1.4  Results 

The  current  work  demonstrates  that  KAM  tori  may  be  constructed  to  represent  orbits 
with  unprecedented  accuracy  given  certain  constraints  and  accurate  enough  trajectory 
information.  Additionally,  largely  separated  formations  of  satellites  with  secular  drifts  on 
the  order  of  0.5  percent  of  the  original  separation  distance  over  10  days  are  demonstrated  for 
multiple  orbit  altitude  and  inclination  combinations  under  the  influence  of  the  geopotential 
expansion  to  order  and  degree  20.  Tight  formations  analyses  yield  secular  drift  rates 
between  satellites  on  the  order  of  4  nanometers  to  1  micrometer  per  second. 
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II.  Background 

This  chapter  will  briefly  investigate  the  current  methods  of  analysis  of  satellite  forma¬ 
tions/constellations.  Then,  exposition  of  the  major  perturbations  experienced  by  satellites 
in  Earth  orbit  will  be  given.  A  description  and  short  history  of  KAM  theory  will  be  pro¬ 
vided.  Finally,  the  chapter  will  end  with  an  introductory  explanation  of  the  formulation 
of  the  dynamics  used  throughout  the  rest  of  the  work. 

2.1  Satellite  Formations  and  Relative  Motion 

One  of  the  first  methods  developed  for  the  description  of  satellite  relative  motion, 
and  that  which  is  predominantly  used  in  the  design  of  rendezvous  missions  and  satellite 
formations,  involves  using  the  Hill-Clohessy-Wiltshire  (HCW)  equations,  created  originally 
for  the  Gemini  Program  [2].  These  equations  are  given  in  the  relative  frame  as 


X  —  2ioy  —  Suj'^x  =  fx 

(2.1) 

y  +  2u}x  =  fy 

(2.2) 

Z  +  =  fz 

(2.3) 

where  x,  y  and  z  are  the  coordinates  in  the  target-centric  frame,  lv  is  the  orbital  angular 
velocity  of  the  chief  satellite  (which  is  equivalent  to  the  mean  motion  n  in  this  circular 
case)  and  the  /„  are  external  forces  [19].  These  equations  are  the  result  of  the  investigation 
of  the  relative  motion  between  a  ‘chief’  (or  ‘target’)  spacecraft  and  a  ‘deputy’  (or  ‘chaser’) 
spacecraft,  after  ignoring  nonlinear  gravity  terms  and  assuming  the  chief  is  in  a  circular 
orbit  [19]. 

In  their  original  form,  the  HCW  equations  provide  an  acceptable  depiction  of  space¬ 
craft  relative  motion  to  first  order  when  the  satellites  in  question  are  relatively  close  to 
one  another.  It  is  commonly  known  that  the  most  potent  perturbative  effect  not  included 
in  these  linearized  HCW  equations  is  due  to  the  terms  in  the  geopotential  expansion  con¬ 
taining  potential  forces  higher  than  zeroth  order  two-body  motion.  There  have  been  a 
variety  of  techniques  used  to  combat  this  source  of  error,  ranging  from  designing  orbit 
constellations  which  are  essentially  invariant  to  higher-order  geopotential  terms  [17]  to 
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modifying  the  HCW  equations  to  include  the  higher-order  terms  to  some  degree,  usually 
only  including  the  coarsest  oblateness  effect,  or  J2  zonal  harmonic  term  (c.f.  §2.2.1  for  a 
more  detailed  explanation  of  the  geopotential  field’s  effects).  Some  missions,  such  as  the 
canceled  USAF  TECHS  AT  21,  were  to  use  formations  that  were  “close  enough”  to  remain 
within  the  linear  regime  of  the  HCW  equations  [12] .  Modified  versions  of  the  HCW  equa¬ 
tions,  including  both  Earth  J2  harmonic  effects  and  nonlinear  differential  gravity  effects, 
have  been  developed  (as  in  [6], [18])  and  applied  to  obtain  better  approximations  over  larger 
separating  distances  and  have  allowed  for  more  efficient  orbital  designs  --  i.e.  orbit  designs 
which  require  fewer  on-orbit  corrections  (and,  thus,  less  fuel)  to  stay  within  certain  relative 
position  constraints.  However,  it  is  certainly  preferable,  if  the  possibility  exists,  to  more 
accurately  represent  the  true  dynamics  of  the  satellites’  respective  motions  around  the 
primary.  This  would  enable  the  ‘capture’  of  all  nuances  of  the  primary’s  influence  which 
would  cause  secular  growth  in  the  separation  between  the  orbiting  bodies. 

2.2  Perturbations  in  Earth- orbit 

Satellites  in  orbit  around  the  Earth  experience  a  number  of  perturbations  which 
tend  to  force  their  dynamics  away  from  those  which  can  be  precisely  described  by  typical 
Keplerian  two-body  motion.  The  types  of  perturbations  which  are  dominant  depend  upon 
the  altitude  of  the  satellite  above  the  Earth’s  surface.  Objects  in  low-Earth  orbit  (LEO  -  up 
to  about  800  km  altitude)  mostly  experience  conservative  perturbative  forces  due  to  Earth’s 
nonsphericity,  along  with  the  non-conservative  force  due  to  atmospheric  drag.  Conversely, 
objects  in  mid-Earth  orbit  (MEO  -  800km  to  30,000km  altitude)  and  geosynchronous  orbit 
(GEO  -  35,780km  altitude)  experience  relatively  little  perturbative  force  from  Earth’s 
nonsphericity,  as  the  Keplerian  term  dominates  at  high  radii,  and  relatively  little  resistive 
force  from  atmospheric  drag,  which  falls  off  exponentially  with  altitude.  Rather,  as  a 
satellite’s  altitude  increases,  its  perturbations  away  from  Keplerian  become  dominated  by 
third-body  effects  from  the  Sun  and  Moon  and  by  the  non-conservative  force  imparted 
by  as  solar  radiation  pressure.  This  thesis  will  focus  solely  on  satellites  in  LEO  and 
mainly  upon  the  effects  of  the  geopotential,  although  atmospheric  drag  will  be  discussed 
in  brief;  hence,  the  dynamics  of  bodies  in  MEO,  GEO  and  above  are  beyond  the  scope 
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of  this  thesis;  however,  the  methodology  developed  and  applied  here  may  be  extended  to 
arbitrary  potential  fields  which  meet  the  KAM  criteria  as  mentioned  below  and  in  Chapter 
III,  and  would  hence,  in  theory,  be  applicable  to  Earth  orbit  classes  including  third  body 
effects,  given  adequate  trajectory  knowledge. 

2.2.1  The  Earth’s  Geopotential.  This  subsection  provides  a  general  overview  of 
the  geopotential  field  through  discussion  of  the  terms  in  the  geopotential  expansion.  For 
a  full  treatment  of  the  geopotential,  c.f  [19]  and  [23] . 

The  derivation  of  the  geopotential  expansion  around  a  solid  body  may  begin  with 
the  familiar  Poisson’s  equation  for  the  gravitational  potential  V : 

VV  =  AttGp  (2.4) 

where  G  is  the  gravitational  constant,  p  is  the  density  of  the  body  and  =  V  •  V  is  the 
standard  Laplacian  operator.  After  evaluating  Poisson’s  equation  in  spherical  coordinates 
(c.f.  [23])  for  positions  outside  of  the  gravitating  body,  the  potential  becomes 
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which  is  termed  the  “expansion  of  the  geopotential  in  spherical  harmonics” ,  and  in  which 
p,  is  the  gravitational  parameter,  r  is  the  radius  of  the  satellite  from  the  Earth’s  center, 
is  the  radius  of  the  Earth,  n  and  m  are  the  degre  and  order  of  the  expansion  (respectively) , 
P™  are  the  associated  Legendre  polynomials,  Cnm  and  Snm  are  the  gravity  field  coefficients 
given  by  an  Earth-gravity  model,  5  is  the  geocentric  latitude  and  A  is  the  east  longitude. 
The  potential  given  in  Eq.  (2.5)  may  be  most  intuitively  considered  as  a  combination  of 
three  types  of  effects:  zonal,  sectoral  and  tesseral  harmonics. 

2. 2. 1.1  Zonal  Harmonies.  Zonal  harmonics  are  encountered  when  order  m 
=  0  and  1  <  n  <  Umax-  The  first  zonal  harmonic,  n  =  1,  has  the  effect  of  moving  the 
center  of  mass  of  the  geoid  north  or  south;  as  such,  it  is  not  usually  included  in  Earth 
geopotential  models,  since  we  desire  the  center  of  our  coordinate  system  to  coincide  with 
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the  center  of  mass.  The  first  nonzero  zonal  harmonic  term,  then,  occurs  when  n  =  2;  this 
term  accounts  for  the  direct  oblateness  (or  equatorial  bulge)  of  the  Earth  due  to  its  spin 
about  the  3-axis.  The  potential  associated  with  n  =  2,  m  =  0  then  assumes  the  form 

V2O  =  ^^|^(3cos20-1),  (2.6) 

where  J2  is  the  C20  coefficient,  with  a  value  ~  0.001082.  This  effect  is  the  second  largest 
in  magnitude,  with  a  value  of  approximately  one  thousandth  of  the  Newtonian  potential 
term  (where  n,m  =  0).  Increasing  the  order  while  maintaining  a  zero  degree  repre¬ 
sents  zonal  harmonic  functions  with  successively  more  “nodes”,  which  in  turn  allows  for 
the  increasingly  refined  modeling  of  mass  distribution  irregularities.  This  idea  of  zonal 
harmonics  is  shown  pictorially  in  Figure  (2.1),  with  both  side  and  top  views,  for  n  =  [2, 6]. 

2. 2. 1.2  Sectoral  Harmonics.  The  sectoral  harmonics  describe  the  effect 
that  occurs  when  n  =  m.  The  associated  Legendre  polynomials  P”  [sin((i)]  have  zeros 
only  when  5  =  ±7r,  i.e.  when  the  satellite  is  over  the  poles.  The  field  terms  cos{m\) 
and  sin{m\)  are  zero  at  2m  lines  of  longitude  A  around  the  geoid;  hence,  the  sectoral 
harmonics  divide  the  globe  into  slices,  much  like  a  typical  beach  ball’s  colored  panels.  The 
n,  m  =  1  term  is  set  to  zero,  since  it  has  the  effect,  as  the  analogous  term  in  the  zonal 
harmonics,  of  shifting  the  planet’s  center  of  mass  away  from  the  center  of  our  coordinate 
frame.  The  (722,  S22  sectoral  harmonic  then  separates  the  earth  into  four  slices:  two  slices 
opposite  each  other  with  higher  potential  and  two  slices  opposite  each  other  with  lower 
potential.  As  was  the  case  regarding  the  zonal  harmonics,  increasing  the  order  and  degree 
divides  the  Earth  into  higher  and  higher  “resolution”  slices.  Figure  (2.2)  shows  sectoral 
harmonics  for  n,m  =  2  and  n,m  =  3. 

2. 2. 1.3  Tesseral  Harmonics.  Tesseral  harmonics  refer  to  the  “off-diagonal,” 
non-zonal  terms:  those  determined  by  n  7^  m,  where  n,  m  >  0.  The  potential  contributions 
from  tesseral  harmonics  are  due  to  mass  distribution  differences  manifested  in  a  tiled  or 
grid  pattern.  In  general,  much  as  before,  the  higher  the  degree  and  order  of  the  expansion. 
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Figure  2.1.  Depiction  of  Zonal  Harmonics  (from  ref.  [19]). 


Figure  2.2.  Depiction  of  Sectoral  Harmonics  (from  ref.  [19]). 


Figure  2.3.  Depiction  of  Tesseral  Harmonics  (from  ref.  [19]). 


the  ‘finer’  the  grid  and  the  higher  the  resolution  of  the  gravitic  model.  Figure  (2.3)  shows 
representative  examples  of  tesseral  harmonics  for  m  =  3,4  and  n  =  [1,3]. 

2.2. 1.4  Earth  Gravity  Model.  To  obtain  a  usable  geopotential  expansion, 
the  Cnm  and  Snm  coefficients  must  be  empirically  determined.  This  has  been  accomplished 
using  satellites  instrumented  with  highly  accurate  equipment  capable  of  detecting  and  doc¬ 
umenting  very  small  perturbations  in  the  satellites’  motion  due  to  differential  gravity  forces. 
Such  a  set  of  coefficients  is  termed  a  gravity  model.  The  current  standard  gravity  model, 
complete  to  order  and  degree  360,  is  called  the  Earth  Gravity  Model  of  1996  (EGM96) 
[19].  The  EGM96  is  highly  accurate  and  has  been  augmented  and  verified  using  data  from 
multiple  satellites.  Although  the  model  includes  terms  up  to  order/degree  360,  the  current 
work  uses  only  terms  to  m,  n  =  20  in  order  to  decrease  computational  burden  while  still 
maintaining  very  high  accuracy;  recall  the  steep  decline  of  the  potential  magnitude  as  order 
and  degree  increase. 

2.2.2  Orbital  Atmospherie  Drag.  As  mentioned  above,  the  main  source  of  dynam¬ 
ical  perturbation  for  satellites  in  LEO  aside  from  the  Earth’s  nonsphericity  is  atmospheric 
drag.  The  current  work  is  concentrated  on  the  KAM  tori  resulting  from  the  conservative 
perturbations  due  to  the  geopotential;  thus,  only  a  brief  discussion  of  the  drag  effects  and 
their  associated  inclusion  into  KAM  theory  will  be  given. 

Atmospheric  drag  manifests  as  a  non-conservative  force  on  a  satellite,  in  a  direction 
opposing  that  of  the  satellite’s  velocity  vector  and  with  a  magnitude  which  depends  on 
several  key  physical  parameters  of  the  satellite  and  its  altitude.  A  standard  form  of  the 
drag  force  equation  on  a  satellite  is  given  as 

fd  =  -lCdAp\Vrel\Vrel  (2.7) 

where  Cd  is  the  drag  coefficient  of  the  satellite,  A  is  its  cross-sectional  area,  p  is  the 
atmospheric  density  at  the  position  of  the  satellite  and  V^ei  is  the  velocity  of  the  satellite 
relative  to  the  local  atmosphere.  Often,  however,  it  is  convenient  to  rewrite  the  drag  effect 
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as  an  acceleration,  rather  than  a  force: 


ad  =  -l^\Vrel\Vrel  (2.8) 

In  Eq.  (2.8),  [3  is  known  as  the  ballistic  coefficient  and  is  given  in  terms  of  the  above 
parameters  and  the  satellite  mass  m  by 


P 


m 


(2.9) 


so  that  a  larger  ballistic  coefficient  is,  of  course,  preferable  when  possible,  since  it  reduces 
the  overall  effect  of  drag  on  the  orbiting  vehicle. 

The  drag  force  on  a  satellite  is  often  one  of  the  most  difficult  aspects  of  the  dynamics 
to  predict,  largely  because  of  the  density  term,  p.  The  density  varies  not  only  with  altitude 
but  can  vacillate  wildly  with  respect  to  solar  activity;  that  is,  in  periods  of  maximum  solar 
flux,  the  atmosphere  reaches  a  state  of  excitement  in  which  the  atmosphere  can  expand 
to  encompass  a  much  larger  volume  and  hence  leads  to  a  greater  apparent  density  at  any 
particular  altitude. 

The  effects  of  drag  on  a  satellite  orbit  can  be  quite  profound.  In  the  case  of  a  satellite 
with  an  elliptical  orbit,  increased  drag  at  its  periapsis  will  cause  a  lowering  of  the  apoapsis 
altitude  in  an  effect  termed  “circularization,”  with  only  a  slight  lowering  of  the  periapsis 
altitude.  This  phenomenon  will  continue  until  the  orbit  is  virtually  circular,  at  which 
time  the  semi-major  axis  a  will  continue  to  decrease,  lowering  the  orbital  altitude  until  the 
satellite  disintegrates  or  impacts  the  planet’s  surface.  Figure  (2.4)  below  shows  the  lifetime 
of  an  example  320km  altitude  orbit  inclined  at  i  =  30°  with  a  nearly  zero  eccentricity,  for  a 
satellite  with  a  ballistic  coefficient  of  /3  ~  318  (corresponding  to  m  =  500kg,  A  ~  0.785m^ 
and  Cd  =  2;  see  §4.3  for  more  detailed  discussion  regarding  inclusion  of  drag  in  numerical 
integration).  The  reason  for  the  nonlinear  decline  after  a  certain  point  in  the  orbit’s  decay  is 
due  to  a  simple  aerodynamic  phenomenon:  as  the  atmospheric  density  grows  exponentially 
with  descent,  the  drag  force  bleeds  energy  from  the  system  and  causes  the  satellite  to  lower 
ever  faster  until  impact  . 
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Figure  2.4.  Example  of  orbital  decay  due  to  atmospheric  drag  for  a  nearly  circular  orbit 
at  320km,  i  =  30° 

Certain  of  the  orbits  examined  in  the  current  work  are  circular  with  initial  altitudes 
on  the  order  of  320km.  The  author  notes  that  satellites  in  such  orbits  will  not,  without 
further  boosting,  have  lifetimes  significant  enough  to  apply  KAM  theorem.  However,  they 
are  utilized  in  the  current  work  because,  with  respect  to  KAM  modeling  in  the  conservative 
perturbations,  they  represent  a  sort  of  “upper  bound”-  i.e.  a  satellite  at  320km  will  be 
experiencing  the  full  effect  of  the  geopotential  perturbations;  thus,  it  stands  to  reason  that 
if  such  a  “worst  case”  can  be  accurately  modeled  using  KAM  theorem  while  ignoring  drag, 
it  should  be  possible  to  extend  the  methodology  to  less  perturbed,  higher  altitude  orbits 
with  no  loss  of  accuracy.  See  Chapters  III  and  IV  for  details  regarding  the  methodology 
and  results  for  varying  orbit  parameters. 

2.3  KAM  Theory 

The  theorem  posited  by  Kolmogorov  [8],  and  later  proved  by  Arnold  [1]  and  Moser 
[15]  possesses  the  possibility  for  far-reaching  application  in  classical  mechanics.  Celletti 
and  Chierchia  have  shown  that  KAM  theory  serves  to  adequately  describe  the  dynamics  of 
celestial  bodies  [3]  and  have  also  applied  it  to  the  dynamics  of  the  Sun-Jupiter-Victoria  sys- 
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tern,  which  is  a  Restricted,  Circular,  Planar  Three  Body  Problem  (RCPTBP)  [4].  McGill 
and  Binney  have  demonstrated  development  of  KAM  tori  to  model  dynamical  systems  in 
generalized  potentials  with  promising  results,  as  in  [13].  Kaasalainen  and  Binney  have 
shown  KAM  torus  development  in  Stackel  potentials  and  in  loop  and  box  orbits  [7].  Ad¬ 
ditionally,  KAM  theory  has  been  utilized  in  the  field  of  quantum  mechanics  to  describe 
the  motion  of  particles  in  magnetic  potentials,  such  as  those  experienced  in  a  particle 
accelerator  (c.f.  [5],  [20], [21]  and  [22]). 

KAM  theory  has  newly  been  applied  to  Earth-orbiting  satellites  by  Wiesel  in  [24]  and 
[25]  and  Little  in  [11].  Wiesel  demonstrated  a  least-squares  method  for  obtaining  KAM  tori 
from  numerically  integrated  data  in  [24],  where  he  showed  the  torus  construction  for  an 
Earth  satellite,  and  later  a  refined  method  using  Eourier  analysis  in  [25],  where  he  showed 
the  construction  of  a  torus  for  a  restricted  three  body  problem  resembling  the  Earth-Moon 
system.  Little  employed  observed  data  from  the  GRACE  and  Jason- 1  satellites  to  show 
that  Earth  satellites  likely  lie  on  KAM  tori  [11]. 

As  mentioned  in  §1.2,  KAM  theory,  at  its  core,  concerns  the  dynamical  behavior  of 
a  system  describable  as  an  integrable  Hamiltonian  subject  to  some  small  perturbation;  i.e. 

H,{I,  ip)  =  Ho{I)  +  if)  (2.10) 

where  is  the  perturbed  Hamiltonian,  I  and  (p  are  the  coordinates  in  an  Action-Angle 
representation,  Hq  is  the  integrable  Hamiltonian,  Hi  is  the  perturbing  Hamiltonian  and  e 
is  some  small,  real  value.  Hq  and  Hi  must  be  smooth,  real-analytic  functions.  According 
to  KAM  theory,  the  motion  of  a  body  in  a  system  which  satisfies  these  conditions  will  lie 
on  a  torus  in  the  system’s  phase  space;  i.e.  the  motion  is  ‘constrained’  to  some  locus  of 
points  in  the  phase  space.  The  dimension  of  the  torus  depends  upon  the  number  of  degrees 
of  freedom  of  the  system:  Hamilton- Jacobi  theorem  mandates  that,  for  a  Hamiltonian  with 
N  degrees  of  freedom,  the  associated  KAM  torus  is  Wdimensional  and  occupies  a  phase 
space  of  2N  dimensions. 
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2.3.1  On  The  Visualization  of  Tori.  One  may  think  of  tori  topologically  as  the 
products  of  circles;  i.e.: 

X  X  ...  X  (2.11) 

' - V - ' 

n 

where  5  is  a  topological  space  which  is  the  subspace  {{xi,X2)  G  :  xf  +  =  l}  of 

Tori  of  dimensions  1  and  2  are  then  rather  easily  visualized;  a  1-torus  is  simply  a  circle 
-  a  1-dimensional  torus  existing  in  2-dimensional  space.  For  example,  is  depicted  in 
Figure  (2.5)  below.  The  2-torus,  then,  is  a  shape  similar  to  a  donut:  a  2-dimensional  torus 
existing  in  4-dimensional  space,  for  one  needs  four  parameters  to  define  a  point  on  its 
surface-  two  actions  and  two  angles  (pk  (c.f.  figure  (2.6)).  KAM  tori  manifest  in  shapes 
topologically  similar  to  this  when  one  examines  systems  with  two  degrees  of  freedom,  such 
as  the  restricted  three-body  problem. 


Figure  2.5.  Example  of  a  1-torus. 

Tori  of  higher  dimensions  (T”,n  >  3)  are  much  more  difficult  for  us  to  visualize, 
and  so  we  must  be  satisfied  with  describing  them  almost  entirely  in  the  language  of  math¬ 
ematics.  However,  one  can  obtain  a  very  basic  idea  of  system  behavior  by  investigating 
higher-dimensional  tori  using  a  method  loosely  reminiscent  of  Henri  Poincare’s  well-known 
“method  of  sections:”  in  the  case  of  an  invariant  3-torus  (such  as  those  examined  in  the 
current  work),  if  one  ignores  one  of  the  angle  coordinates  and  its  associated  action,  one 
can  obtain  a  two-dimensional  “projection”  of  the  3-torus  in  the  form  of  2-torus.  Of  course. 
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Figure  2.6.  Top(left)  and  perspective  (right)  views  of  a  2-torus. 

the  amount  of  information  one  gains  from  such  an  exercise  is  only  useful  insomuch  as  it 
provides  a  very  general  idea  of  torus  action  proportions  and  frequency  trends.  Shown 
below  in  Figure  2.7  is  such  a  2-torus  “projection”  obtained  during  the  present  work  for  a 
KAM  3-torus  of  an  Earth-orbiting  satellite,  where  the  smallest  frequency  and  its  associated 
action  are  ignored. 


Figure  2.7.  Example  of  a  2-toroidal  “projection”  of  the  KAM  3-torus  for  an  Earth¬ 
orbiting  satellite. 
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^•4  Dynamics  Formulation 

In  order  to  choose  an  appropriate  coordinate  frame  in  which  to  formulate  the  system 
dynamics  and  construct  KAM  tori,  one  must  strive  to  find  the  frame  which  affords  the 
largest  number  of  integrals  of  the  motion.  In  the  general  case  (posed  in  Cartesian  coordi¬ 
nates  in  an  inertial  frame)  of  a  satellite  moving  in  Earth’s  full  potential  field,  there  are  no 
integrals  of  the  motion  readily  available.  However,  we  may  find  a  constant  of  the  motion 
in  this  case  by  choosing  a  frame  which  rotates  with  the  Earth-  namely,  the  Earth-centered 
rotating  frame,  also  called  the  Earth-centered,  Earth- fixed  (ECEE)  frame,  in  which  the 
1-axis  points  from  the  center  of  the  earth  through  the  intersection  of  prime  meridian  and 
the  equator,  the  3-axis  points  along  the  axis  of  rotation,  and  the  2-axis  completes  the 
right-handed  orthonormal  basis.  In  the  ECEE  frame,  the  elements  of  nonsphericity  which 
cause  the  differential  gravitational  contributions  in  the  Earth’s  potential  are  fixed;  i.e.,  the 
position  vector  of  a  specific  infinitesimal  piece  of  the  Earth’s  mass  will  always  remain  con¬ 
stant.  We  take  advantage  of  this  fortunate  fact  to  greatly  simplify  our  equations  of  motion. 
We  begin  our  Hamiltonian  formulation,  following  Wiesel  [24],  by  writing  the  expressions 
for  the  specific  momenta  pn  as 


Px  =  x-iV(£,y  (2.12) 

Py  =  y  +  U(£,x  (2.13) 

Pz  =  z  (2.14) 

where  x,  y  and  z  are  the  coordinates  of  the  satellite  in  the  ECEE  frame  and  Wq  is  the 
angular  rate  of  rotation  of  the  Earth.  We  may  then  find  the  Lagrangian  for  the  system 

L  =  T-V  (2.15) 

where  T  and  V  are  the  kinetic  and  potential  energies,  respectively.  We  can  use  the  La¬ 
grangian  to  calculate  the  Hamiltonian  using: 


-  L 

i 


(2.16) 
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After  substituting  the  appropriate  expressions  for  kinetic  and  potential  energies,  the  Hamil¬ 
tonian  may  be  written  in  its  final  form 


H  =^{pI+pI+  pI)  +  uJ(B{yPx  -  xpy) 

oo  n 

n=0  m=0  ® 

X  [CnmCOs{mX)  SnmSin{mX)] 


(2.17) 


where  p  is  the  gravitational  parameter,  r  is  the  radius  of  the  satellite  from  the  Earth’s 
center,  i?0  is  the  radius  of  the  Earth,  are  the  associated  Legendre  polynomials  and  Cnm 
and  Snm  are  the  gravity  field  coefficients  given  by  an  Earth-gravity  model,  as  described  in 
the  previous  section.  The  radius  r,  the  geocentric  latitude  5  and  the  east  longitude  A  of 
the  satellite  may  be  found  by  elementary  geometry  to  be: 


sin6 

tanX 


z 

^3,2  _|_  y2 

y 

X 


(2.18) 


After  following  the  above  method,  we  may  obtain  the  one  achievable  constant  of  the  motion 
in  this  problem:  the  Hamiltonian.  Again  we  note  that  this  approach  includes  only  the 
potential  perturbations  and  ignores  all  non-conservative  forces. 

Throughout  this  work,  a  variety  of  units  and  constants  are  used  in  calculations  and 
resulting  displays.  Most  of  the  units  are  common,  such  as  the  standard  degree,  radian, 
meter  and  kilometer.  A  few  are  not  quite  so  usual,  such  as  the  time  unit  (TU)  and  the 
Earth  radius  (ER).  The  values  of  these  units  and  of  assorted  constants  used  throughout 
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this  work  are  listed  below  [19].  Minutes  and  seconds  given  are  in  solar  time. 

1  TU  =  13.446852  min 
1  ER  =  6378.137  km 

w©  =  0.0588335998  rad/TU 

nlTl  ^ 

//©  =  398600.4418^^  =  1  ER^/TU'^ 
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III.  Method 


The  core  of  utilizing  KAM  theory  for  orbital  mechanics  lies  in  calculating  the  KAM  tori  for 
orbits  in  question.  This  chapter  will  begin  with  a  discussion  of  the  process  of  application 
of  KAM  theorem  to  satellite  orbits.  It  will  then  proceed  with  details  of  the  methods  used 
to  obtain  the  requisite  propagated  orbital  data.  Finally,  an  explication  of  methodology 
used  for  formation  analysis  on  the  torus  will  take  place. 

3.1  Describing  Orbits  as  KAM  Tori 

As  mentioned  before,  the  KAM  torus  for  an  Earth  satellite  exists  as  a  3-torus  in 
6-space  in  terms  of  three  angle  variables  and  their  associated  actions.  However,  in  order 
to  use  orbital  data  to  gain  any  utility  from  the  concept  of  a  KAM  torus,  we  must  find  a 
formulation  of  the  torus  which  involves  the  physical  coordinates.  Ideally,  what  we  would 
seek  is  a  map 

M  :  {q,p)  ^  (3.1) 

which  would  allow  us  ultimately  to  find  the  torus  actions  Ik  and  angles  (pk^  having  as  the 
input  a  state  containing  the  physical  coordinates  qk  and  momenta  Unfortunately,  there 
does  not  seem  to  exist  so  simple  a  homeomorphism  between  the  physical  space  and  torus 
space.  To  circumvent  this  problem,  we  must  obtain  the  states  at  a  number  of  points  along 
the  orbital  path  (through  satellite  observation  or  numerical  propagation)  and  perform  an 
analysis  upon  these  data  to  find  two  separate  types  of  information:  the  fundamental  or 
basis  frequencies  with  which  the  system  oscillates  (which  are  related  to  the  torus  angle 
variables  along  with  their  temporal  development)  and  the  amplitudes  with  which  each  of 
the  frequencies  and  combinations  thereof  manifest  in  the  orbital  motion  (related  to  the 
torus  actions).  In  an  isoenergetic  system  (one  where  the  only  forces  are  conservative,  or 
potential,  forces),  the  torus  actions  I  will  be  constant  -  indeed,  one  of  the  beauties  of 
representing  the  dynamics  as  a  torus  is  that  the  Ik  are  constant  and  the  pk  increment 
linearly  with  time.  What  is  found  through  the  below  method,  in  fact,  is  the  reverse  of  the 
map  described  above,  manifested  as  a  series  involving  torus  parameters  used  to  calculate 
physical  coordinates  and  then  their  associated  momenta. 
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3.1.1  KAM  Series  Representation.  According  to  KAM  theorem  and  Hamilton- 
Jacobi  theorem,  a  nearly  integrable  system  of  3  degrees  of  freedom  (such  as  that  in  the 
current  problem)  should  also  have  three  fundamental  frequencies,  which,  it  has  been  es¬ 
tablished,  are  intimately  related  to  the  representative  angle  coordinates  of  the  tori.  Due  to 
the  conservative  and  nearly  integrable  nature  of  this  system,  and  since  we  seek  a  function 
that  is  periodic  in  the  angle  variables  99^,  it  is  natural  to  seek  a  representation  of  the  orbit 
coordinates  as  a  Fourier  series: 


(3.2) 

j 

Again,  what  we  work  to  find  through  frequency  analysis  will  be  the  complex  series  coeffi¬ 
cients  Dj  and  their  associated  frequency  combinations  j  -(p.  The  ipk  can  be  described  using 
the  basis  frequencies  and  time  t  through 

tpk  —  -|-  (3-3) 

where  is  the  /cth  fundamental  frequency  and  p>kQ  is  the  initial  fcth  angle  value  {(pk  at 
t  =  0).  In  Eqn.  3.2,  each  j  is  a  vector  of  integers,  which,  when  the  inner  product  j  ■  ip 
is  taken,  allows  for  specification  of  integer  combinations  of  the  different  basis  frequency 
contributions.  The  details  and  mechanics  of  this  will  be  discussed  below. 

3.1.2  Speetral  Analysis.  The  problem  of  finding  the  basis  frequencies  for  an  orbit 
is  best  solved  using  a  method  similar  to  that  of  Laskar  as  given  in  [9]  and  [10],  known 
as  the  Numerical  Algorithm  of  the  Fundamental  Frequency  (NAFF).  Per  Laskar,  we  may 
take  the  finite  fourier  transform  of  coordinate  data  f{t)  at  a  frequency  uj  by 

=  ^  fj{t)e-^^\{t,T)dt  (3.4) 

Laskar’s  method  is  obviously  quite  similar  to  normal  Fourier  analysis,  but  uses  instead  a 
weight  or  window  function  x{t)  to  circumvent  accuracy  losses  caused  by  the  fact  that  the 
orbital  data  f{t)  is  likely  not  truncated  perfectly  to  contain  an  integer  number  of  periods 
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over  the  time  span  [— T,  T] .  The  aforementioned  losses  are  due  to  a  phenomenon  known 
as  frequency  leakage-  c.f.  §3. 1.2. 2. 

3. 1.2.1  The  Power  Spectral  Density.  The  goal  of  Laskar’s  method,  when 
searching  for  a  basis  frequency,  is  to  find  the  maximum  of  Eqn  (3.4)  in  the  neighborhood 
of  a  rough  initial  guess  of  the  basis  frequency  in  question.  This,  in  essence,  allows  for  the 
determination  of  the  frequency  value  (accurate  to  within  computational  uncertainties)  at 
which  the  spectral  power  is  the  highest.  The  spectral  power,  defined  as  V  =  is  the 
most  common  way  to  quantify  contributions  from  different  frequencies  inherent  in  a  signal. 
A  plot  of  the  power  values  for  a  range  of  frequencies,  or  spectrum,  is  aptly  termed  a  Power 
Spectral  Density  (PSD)  plot  and  is  a  useful  way  of  rather  intuitively  examining  the  periodic 
characteristics  of  the  signal.  An  example  of  a  PSD  plot  for  an  orbit  is  shown  in  Figure 
(3.1)  below.  By  simple  observation,  it  is  obvious  for  a  non-chaotic  orbit  that  there  are  a 
number  of  peaks  in  the  PSD  at  various  frequencies.  These  power  peaks  correspond,  in  the 
case  of  our  perturbed  orbit  data  “signal”,  to  combinations  of  the  three  basis  frequencies. 


0  0.5  1  1.5  2  2.5  3 

Frequency  (rad/s) 


Figure  3.1.  Example  of  Power  Spectral  Density  plot  of  an  orbit  for  w  =  [0,7r]  with 
window  power  p  =  2. 
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3. 1.2.2  Window  Function  and  Window  Power.  The  window  function  used 
in  the  current  work  for  T)  in  Eqn.  (3.4)  is  known  as  the  Hanning  window,  and  is  given 
by 


x{t,T)  = 


2P  {p'.y 
(2p)! 


1  f 

1  +  cos  (  — 


(3.5) 


where  p  is  the  window  power  and  t  and  T  are  the  current  time  step  and  the  end  time, 
respectively,  as  discussed  above.  Essentially,  this  window  function  serves  to  “smooth”  the 
ends  of  the  data  set  so  that  the  values  taper  down  to  zero  with  varying  speeds  (related  to 
the  values  of  p) .  Such  a  tapering  of  the  data  eliminates  the  frequency  leakage  mentioned 
above  by  eliminating  the  discontinuity  encountered  if  one  “ties”  the  beginning  and  end  of 
the  data  set  together  and  treats  it  as  one  continuous  periodic  system.  Laskar  demonstrates 
efficient  use  of  the  Hanning  window  with  powers  of  3  — 5  and  states  that  powers  above  p  =  5 
tend  to  decrease  accuracy  [9];  Wiesel  [24]  and  Little  [11]  use  p  =  2  almost  exclusively.  The 
author  has  found  in  the  current  work  that  a  general  statement  about  appropriate  values 
for  p  cannot  be  made;  that  is,  the  choice  of  power  depends  upon  the  accuracy  of  the  signal 
data,  the  length  of  the  signal  sampling  period  and  the  “stiffness”  of  the  dynamics  (related, 
in  this  case,  to  the  order  and  degree  of  the  geopotential  expansion  considered),  among 
other  factors.  Wiesel  shows  an  instructive  example  of  the  effects  on  the  PSD  of  a  spectral 
line  by  increasing  window  powers  in  [25] ;  essentially,  as  one  increases  the  value  of  p,  the 
values  of  the  sidelobes  around  the  spectral  line  fall  off  more  and  more  quickly,  with  the 
potential  sacrifice  of  frequency  accuracy  due  to  the  broadening  of  the  true  frequency’s  peak 
in  the  PSD.  See,  for  example,  Eigure  (3.2)  below.  As  in  Wiesel,  a  value  of  p  =  2  is  a  good 
general-purpose  power;  however,  as  mentioned  by  Laskar  [9] ,  it  is  advantageous  to  explore 
higher  window  powers  until  accuracy  gains  cease.  In  the  current  work,  a  power  p  =  2 
proved  sufficient  for  analysis  of  orbit  integrations  including  n,m  <  5;  however,  for  1-year 
integrations  with  n,m  =  20,  gains  in  accuracy  of  the  fundamental  frequencies  continued 
in  some  cases  up  to  p  =  7,  allowing  the  to  be  determined  to  within  an  error  of  a  few 
parts  in  10^^. 


3. 1.2.3  Fourier  Coefficients.  As  previously  mentioned,  the  torus  actions  Ik 
are  indirectly  reflected  in  the  complex  Eourier  coefficients  Dj  in  Eqn.  (3.2).  It  is  sometimes 
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Figure  3.2.  Effects  on  the  PSD  for  a  spectral  line  of  a;  =  0.057  rad  from  increasing 
Hanning  window  power  p  in  Eqn.(3.5)-  based  on  example  from  Wiesel,  [25] 


advantageous  to  rewrite  the  series  (3.2)  in  real  form  as 


q{t)  =  '^  [Cjcos  (j  •  if)  +  S-jSin  {j  ■  p)]  (3.6) 

j 

so  that  we  seek,  from  our  spectral  analysis,  the  coefficients  and  obtainable  from 
Eqn.  (3.4)  by 


Ck  =  23f?<(.  {taj)  (3.7) 

Sk  =  (a;  •)  (3.8) 

where  cuj  is  the  frequency  combination  (or  the  PSD  bin  associated  with  the  j  in  question), 
calculated  with 

io-j=j-n  (3.9) 

The  constant  term  for  each  coordinate,  Cq^,,  may  be  found  by  evaluating  Eqns.  (3.7)  and 
(3.8)  with  a;  =  0. 
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3. 1.2. 4  Spectral  Decomposition.  Finding  the  values  of  the  Fourier  coef¬ 
ficients  mentioned  in  §3. 1.2. 3  has  proven  to  be  somewhat  of  an  art.  The  author  of  the 
current  work  has  found  that  simple  Fourier  analysis,  i.e.  obtaining  the  amplitudes  of  the 
spectral  peaks  per  Eqns.  (3.7)  and  (3.8),  is  only  effective  to  a  certain  accuracy  level  due  to 
two  phenomena:  spectral  “shadowing”  and  frequency  folding.  Spectral  shadowing  occurs 
when  one  is  dealing  with  a  system  where  there  is  one  very  small  frequency  among  the  basis 
frequency  set;  i.e. 

3  {Dm,  Dn)  G  D  I  Dm  <C  Dn,  m  /  n  (3.10) 

where  the  frequency  proportion  considered  “small”  depends  upon  the  accuracy  of  the  signal 
data  and  computational  precision.  The  effect  this  has  on  the  spectral  analysis  can  be  quite 
profound,  because  when  one  is  analyzing  the  amplitudes  of  integer  combinations  of  the 
elements  of  D,  there  will  be  circumstances  in  which  one  is  interested  in  some  frequency 
combination  involving  the  largest  frequency  plus  or  minus  the  smallest  frequency;  for 
example,  in  a  [not  uncommon]  case  where  we  have  an  index  and  basis  frequency  set  such 
that 


j  =  [101]  (3.11) 

D  =  [0.98  0.060  0.0018]  rad/s  (3.12) 

we  would  then  then  be  seeking  the  amplitude  of  a  frequency 

COj  =  j  •  D 

=  0.9818 

In  this  case,  the  power  of  the  basis  frequency  at  fli  =  0.98  will  almost  certainly  have 
a  higher  amplitude  than  coj,  and  since  loj  is  so  near  to  fli,  it  will  be  very  difficult  to 
distinguish  the  true  value  of  the  amplitude  at  tOj.  This  issue  is  even  further  exacerbated 
when  using  a  window  function,  which,  as  mentioned  above  in  §3. 1.2. 2  and  shown  in  Fig. 

(3.2),  has  the  effect  of  widening  the  peak,  thus  further  obscuring  the  contributions  of  any 

frequencies  in  near  proximity. 
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The  second  malicious  effect,  and  an  issue  inherent  in  any  practical  signal  analysis 
application,  is  frequency  folding,  which  involves  the  “reflection”  of  spectral  peaks  about 
the  Nyquist  frequency  or  Nyquist  limit, 


^Nyquist  —  ,  (3.13) 

^sample 

where  tgampie  is  the  sampling  time  interval  of  the  [orbital]  signal  data.  For  example, 
consider  the  simple  periodic  function 


X(t)  =  sin(r'ot)  (3-14) 

where  the  frequency  of  oscillation  i/q  =  0.55  rad/TU  and  we  obtain  “propagation”  data  for 
times  t  =  [— T  :  5t  -.T],  where  we  take  T  =  8000  TU  and  6t  =  0.5  TU.  The  PSD  for  this 
case  after  Fourier  transform  per  Eqn.  (3.4)  is  shown  in  Fig.  (3.3)  over  the  Fourier  domain 
(j  =  [0,  dvr];  it  is  clear  that  the  frequency  domain  representation  of  this  function  is  merely  a 
single  spectral  line  at  cu  =  i^o-  The  reader  is  invited  to  notice,  however,  the  aforementioned 
frequency  folding  manifested  here:  the  signal  in  question  has  a  Nyquist  frequency  of 

TT 

^Nyquist  —  T  —  ^TT  (3.15) 

^sample 

and  it  is  easily  seen  that  there  are  reflections  evident  in  the  PSD  plot  over  the  given 
domain.  One  “false”  peak  occurs  at  ujn  =  oj Nyquist  —  one  at  ^^2  =  ^Nyquist  +  t'o  and 
one  at  ujr^  =  Nyquist  —  The  pattern  of  false  peaks  thus  established  continues  on 
ad  infinitum  as  one  examines  the  powers  of  higher  and  higher  frequencies.  In  general, 
however,  we  are  only  interested  in  those  frequencies  at  or  below  oj  Nyquist-,  as  the  Nyquist 
frequency  represents  the  highest  frequency  that  can  be  accurately  sampled.  We  are  left, 
then,  with  the  issue  of  the  reflection  peak  at  =  oj  Nyquist  —  t'O)  since  it  poses  a  problem 
when  we  are  trying  to  determine  powers  of  the  real  frequency  combinations  inherent  in  the 
dynamical  system. 

Conveniently,  it  seems  that  the  solution  to  both  the  spectral  shadowing  and  fre¬ 
quency  folding  problems  lies  in  a  spectral  decomposition  method.  Spectral  decomposition 
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Figure  3.3.  PSD  of  a  spectral  line  at  i^o  =  0.55  rad/TU  and  sampling  time  6t  =  0.5  TU 
demonstrating  frequency  folding. 

consists  of  proceeding  through  the  spectral  content  of  a  signal,  starting  with  the  frequency 
combination  having  the  largest  power  V{ui),  and  removing  the  signal  contribution  of  each 
spectral  line  before  proceeding  to  the  next  highest.  In  essence,  if  one  considers  the  fre¬ 
quencies  and  their  associated  powers  composing  the  spectrum  of  a  signal,  organized  into  a 
set, 

W  =  [UJI,UI2,  .  .  .,UJm]  I  V{iOn+l)  <  V{uJn)  (3.16) 

where  m  is  the  total  number  of  frequencies  to  be  analyzed  in  the  spectrum,  the  spectral 
decomposition  process  proceeds  iteratively  through  W  to  find  the  complex  amplitudes 
using  the  equation 

4>n  =  4>{^n)  =  ^  J  (3.17) 

where  fn{t)  represents  the  “decomposed”  orbital  signal  data,  having  removed  from  the 
original  data  /o  the  contributions  of  all  frequencies  w  <  in  W;  that  is, 

n—1 

(3.18) 

h=l 
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In  the  case  of  spectral  shadowing,  the  use  of  this  decomposition  method  allows  the  analyst 
to  distinguish  the  separate  contributions  from  two  proximate  frequencies  with  a  much 
higher  accuracy  than  is  possible  with  mere  Fourier  analysis.  In  the  case  of  frequency 
folding,  spectral  decomposition  may  be  applied  with  one  important  caveat:  in  all  cases 
observed  so  far  by  this  author,  the  “true”  peak  will  have  a  higher  power  in  the  PSD  than 
any  of  its  dependent  reflections  under  lONyquisu  that  is, 

^  ^  ^Nyquist  (3.19) 

This  is  evident  in  Fig.  (3.3)  when  comparing  the  peak  at  uq  to  the  peak  at  =  uj^yquist  — 
uq.  This  relationship  betwen  the  powers  of  a  frequency  and  its  reflections  allows  the  analyst 
to  use  the  iterative  process  given  in  Eqns.  (3.17)  and  (3.18),  since  the  conditions  on  W  as 
given  in  Eqn.  (3.16)  are  met.  In  essence,  after  one  removes  the  contribution  of  the  true 
peak  at  cu  =  z/Q)  the  reflection  peaks  disappear  from  the  remaining  signal  and  cause  no 
further  confusion  in  determination  of  peaks  and  their  frequencies. 

3. 1.2. 5  Fourier  Indiees.  The  spectral  analysis  process  outlined  in  the  above 
sections  requires,  upon  mechanization  in  a  numerical  routine,  a  somewhat  unusual  consid¬ 
eration  regarding  the  integer  indices  j  =  ji  +  j2  +  is-  Following  Wiesel  [25],  the  current 
work  utilizes  an  index  coding  strategy  in  which  the  first  non-zero  term  in  j  is  never  nega¬ 
tive;  this  precaution  serves  to  eliminate  repetition  due  to  the  nature  of  negative  arguments 
in  trigonometric  functions-  namely,  sin{—x)  =  —sin{x)  and  cos{—x)  =  cos{x). 

3.1.3  Nature  of  the  Basis  Frequeneies.  The  fundamental  frequencies  manifest 
as  combinations  of  the  frequencies  of  commonly  observed  effects  in  Earth  orbit  [24].  The 
first  and  largest  fundamental  frequency  is  known  as  the  anomalistie  frequeney  and  is  given 
approximately  by 

where  a  is  the  semi-major  axis,  e  is  the  orbital  eccentricity,  i  is  the  orbital  inclination  and 
J2,  i?©  and  ^  are  as  given  previously  in  Chapter  II,  and  where  only  effects  of  the  first 
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zonal  harmonic  are  included  (hence  its  approximate  nature).  It  is  clear  that  this  is  the 
largest  frequency  due  to  its  containing  the  Keplerian  (or  two-body)  frequency  (also  the 
mean  motion  in  our  circular  case), 


ujtb  = 


(3.21) 


as  its  first  term.  The  second  fundamental  frequency  is  the  combination  of  the  earth’s 
rotational  frequency  with  the  orbit’s  nodal  regression  rate: 

3^J2R%  ,  ,  ,  , 

"^2  ~  +  2a7/2(i 

The  nodal  regression  rate  in  this  equation  refers  to  the  tendency  of  the  orbital  plane 
to  rotate  about  the  primary’s  Earth-centered  inertial  (ECI)  frame’s  3-axis  due  to  the 
oblateness  of  the  primary,  most  easily  seen  as  a  secular  change  in  the  angle  between  the 
inertial  1-axis  and  the  line  of  nodes.  The  third,  final  and  smallest  fundamental  frequency 
is  found  to  be  the  apsidal  regression  rate: 

which  accounts  for  the  tendency  of  the  orbit  to  rotate  within  its  own  plane  about  the  orbit 
normal;  this  manifests  as  a  secular  change  in  the  argument  of  perigee  used  in  the  classical 
orbital  elements. 

Depicted  in  Figures  (3.4)-(3.6)  below  are  the  relationships  between  the  fundamental 
frequencies  (as  approximated  by  Eqns.  (3.20)-(3.23))  and  the  inclination  i  and  altitude 
for  nearly  circular  orbits  (e  ~  0.00108).  The  first  basis  frequency  Di,  shown  in  Fig.  (3.4), 
shows  the  expected  behavior  in  that  the  anomalistic  frequency  should  decrease  in  value 
as  the  orbit  increases  in  altitude,  which  is  a  direct  corollary  of  Eqn.  (3.21)  above.  Addi¬ 
tionally,  for  a  given  altitude,  the  first  basis  frequency  decreases  slightly  as  the  inclination 
of  the  orbit  increases.  This  follows  from  Eqn.  (3.20),  as  the  total  contribution  from  the 
second  term  depends  on  sin{i),  and  the  magnitude  is  thus  fully  subtracted  when  i  =  90°. 
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Figure  3.4.  Behavior  of  Oi  for  a  near-circular  orbit  with  varying  altitudes  and  inclina¬ 
tions. 

Figure  (3.5)  shows  that  ^2  decreases  in  value  for  increasing  altitude;  again,  this 
follows  logically  from  the  fact  that  the  difference  between  the  true  geopotential  and  the 
two-body  potential  decreases  as  radius  increases.  The  figure  also  demonstrates  the  well- 
known  behavior  of  the  orbits’  nodal  regression  rates;  these  rates  decrease  nonlinearly  as 
the  orbit  inclinations  increase,  so  that  the  nodal  regression  is  a  maximum  at  i  =  0°  and 
is  zero  sX  i  =  90°,  in  agreement  with  Eqn.  (3.22).  It  is  important  to  note  that  the 
frequency  shown  along  the  vertical  axis  of  Fig.  (3.5)  is,  as  stated  above,  the  combination 
of  the  earth’s  rotation  rate  ujq,  =  0.0588336rad/TC/  with  the  nodal  regression;  this  explains 
why  the  fundamental  frequency  tends  towards  as  i  — >  90°.  Also  of  note  is  that  the 
nodal  behavior  for  retrograde  orbits  (orbits  with  i  >  90°)  is  not  shown  in  Fig.  (3.5).  For 
retrograde  orbits,  the  nodal  regression  rate  becomes  negative,  which  means  that  the  orbit 
processes  westward  rather  than  eastward  as  in  prograde  orbits. 

The  third  basis  frequency,  fls,  depicted  in  Fig.  (3.6),  shows  another  commonly 
characterized  and  utilized  orbital  trend.  This  apsidal  regression  rate  has  a  so-called  critical 
inclination  at  i*  ~  63.4°;  that  is,  the  argument  of  perigee  essentially  does  not  grow  secularly 
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Figure  3.5.  Behavior  of  ^2  for  a  near-circular  orbit  with  varying  altitudes  and  inclina¬ 
tions. 

for  an  orbit  at  this  inclination.  A  common  historical  usage  of  the  critical  inclination 
has  been  the  Molniya  orbit  used  by  the  former  USSR  and  now  Russia  in  order  to  allow 
maximum  coverage  of  the  northern  hemisphere  by  a  satellite  with  a  highly  elliptic  orbit  (c.f. 
[19]  for  more  details).  Figure  (3.6)  also  shows  the  relationship  of  the  apsidal  regression  rate 
of  circular  orbits  to  varying  altitude  and  inclination  when  i  ^  i* .  When  i  <  i* ,  increasing 
altitude  decreases  the  regression  frequency,  whereas  the  opposite  is  true  for  i  >  i*. 

The  fundamental  frequencies  described  in  the  equations  and  plots  of  this  section 
represent,  again,  merely  approximations  of  the  true  frequencies  of  the  system,  as  they 
include  only  the  J2  term  of  the  geopotential  expansion.  In  order  to  find  the  true  values 
of  the  basis  frequencies,  frequency  analysis  must  be  utilized  on  the  propagated  orbital 
data  including  the  other  geopotential  terms  using  the  methodology  outlined  earlier  in  this 
chapter. 


29 


Figure  3.6.  Behavior  of  Os  for  a  near-circular  orbit  with  varying  altitudes  and  inclina¬ 
tions. 

3. 2  Orbital  Propagation 

The  principles  of  §3.1  make  it  clear  that  the  beginning  of  KAM  torus  construction 
lies  in  the  acquisition  of  orbital  data  in  some  fashion.  Judging  from  the  literature,  there  has 
been  considerable  success  in  this  area.  As  mentioned  in  Chapter  II,  Little  has  demonstrated 
fitting  of  KAM  tori  to  observed  orbital  data  for  Earth  satellites  with  reasonable  accuracy 
[11].  Additionally,  Wiesel  has  demonstrated  incontrovertible  evidence  that  a  general  Earth 
satellite  lies  on  a  KAM  torus  through  the  least-squares  fitting  of  numerically  integrated 
data  [24]  and  also  shown  the  fitting  of  numerical  data  for  the  restricted  three-body  problem, 
a  system  with  2  degrees  of  freedom,  using  an  improved  method  of  Fourier  Analysis  [25]. 
The  current  work  builds  upon  the  methodology  used  by  Wiesel  and  Little  and  extends  it 
to  formation  design. 

3.2.1  Equations  of  Motion.  Orbital  data  for  the  orbit  examined  in  the  current 
work  are  obtained  through  numerical  integration  of  the  equations  of  motion  derived  from 
the  Hamiltonian  derived  in  Chapter  II  with  orbit  geometry  determined  by  judicious  choice 
of  initial  conditions.  The  equations  of  motion  may  easily  be  derived  from  the  Hamiltonian 
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by: 


dH 

=  tt- 
OPk 

f)  M 

Pk  =  -^  (3.24) 

OQk 

The  derivative  of  the  state  vector  at  any  position  then  becomes,  after  substitution  and 
differentiation, 


C{q,p)  =  [q  pf 


Pi  +  ^®q2 
P2  - 


P3 


W©P2  — 


UjQsPl  — 


dUjq) 

dqi 

dUjq) 

dq2 


dU{q) 

dq3 


(3.25) 


where  U{q)  represents  the  potential  function  at  position  q  (as  described  in  Ch.  II),  which 
implies  that  its  derivative  represents  the  geopotential  force  upon  the  satellite. 


3.2.2  Numerical  Integrator.  The  numerical  integration  for  this  work  was  per¬ 
formed  using  an  explicit  8th  order  Runge-Kutta  integrator,  based  upon  that  given  by 
Dorman  and  Prince  in  [16],  which  performs  13  evaluations  per  time  step  and  has  a  local 
error  of  order  .  This  integrator,  like  the  Hamming  integrator  used  by  Wiesel  in  [24]  and 
by  Little  in  [11],  is  not  symplectic;  that  is,  it  does  not  explicitly  conserve  the  Hamiltonian 
(this  system’s  constant  of  the  motion).  This  fact  allows  for  us  to  ‘check’  the  accuracy  of 
the  integrator  by  calculating  at  every  output  time  step  the  error  value 


AH  =  H{t)  -  H{to)  (3.26) 

and  ensuring  that  this  difference  has  remained  suitably  small.  Figure  (3.7)  below  shows 
the  error  growth  for  an  archetypal  numerical  integration  of  an  orbit  at  320  km  altitude, 
inclined  at  30  degrees,  over  a  total  time  period  of  approximately  1  year  and  including  all 
geopotential  terms  up  to  order  and  degree  m,n  =  20.  As  seen  in  the  figure,  the  error 
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for  this  typical  case  shows  pseudorandom  behavior  with  maximum  amplitude  near  the 
accuracy  bounds  of  machine  double  precision  at  approximately  4  x  10“^^,  establishing 
a  high  confidence  in  the  accuracy  of  the  integration  results.  Figure  (3.8)  shows  the  same 
Hamiltonian  error  array  after  the  application  of  the  Hanning  window  (Eqn.  (3.5))  with  p  = 
6.  This  represents  the  apparent  error  seen  by  the  spectral  analysis  algorithm.  Comparison 
of  the  two  figures  reveals  two  important  realizations:  first,  the  maximum  error  is  decreased 
from  approximately  4  x  10“^^  to  approximately  2.4  x  10“^^;  second,  the  tradeoff  comes  in 
that  the  local  error  closer  to  the  center  of  the  windowed  data  set  is  slightly  higher  than 
the  associated  error  in  the  unwindowed  data. 


Figure  3.7.  Example  of  Hamiltonian  error  over  integration  time  of  approx.  1  year  for 
an  Earth  orbit  with  m,n  =  20,  320km  altitude,  30  deg  inclination. 

3.2.3  Integration  Characteristics.  The  method  of  spectral  analysis  (and  thereby 
KAM  torus  construction)  outlined  above  requires  integrated  data  with  certain  characteris¬ 
tics  to  function  properly  and  with  the  desired  accuracy.  Eirst,  the  integration  time  must  be 
symmetrically  split  around  t  =  0;  that  is,  the  final  set  of  data  consists  of  a  backwards  inte¬ 
gration  over  the  timescale  — T  <t<0  combined  with  a  forward  integration  over  0  <  t  <  T, 
consistent  with  Eqn.  (3.4).  This  formulation,  as  described  in  [9]  and  [25]  allows  the  user  to 
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Figure  3.8.  Hamiltonian  error  as  in  Fig.  (3.7)  after  application  of  a  window  function 
with  p  =  Q,  per  Eqn.  (3.5). 

take  advantage  of  the  finite  accuracy  available  from  any  real-world  numerical  integration, 
as  the  total  error  growth  over  the  interval  [0,  T]  is  (in  these  and  similar  systems)  smaller 
than  the  error  growth  over  the  interval  [0,2r].  Second,  and  unfortunately,  the  method 
of  KAM  torus  construction  as  it  currently  stands  requires  integration  of  long  periods  of 
orbital  data  (large  T)  with  very  fine  temporal  resolution  in  the  outputs  (small  6t). 

3.3  Formation  Analysis  on  the  Torus 

Using  KAM  theorem  to  design  a  satellite  formation  or  constellation  requires,  of 
course,  that  the  designer  obtain  the  KAM  torus  for  a  particular  orbit  by  applying  the 
methodology  previously  outlined  in  this  chapter  (§3.1)  to  orbital  data,  either  observed  or 
integrated  (per  §3.2).  After  the  torus  series  given  in  Eqns.  (3.2)  and  (3.6)  is  known  in  terms 
of  the  series  coefficients  and  and  the  basis  frequencies  the  astrodynamicist  can 
displace  the  initial  condition  of  a  satellite  while  constraining  it  to  the  same  torus  by  using 
the  equation: 

ipl  =  (pf.  +  Spk  (3.27) 


33 


where  (pk  is  the  kth  angle  variable  as  defined  in  Eqn.  (3.3).  In  effect,  Eqn.  (3.27)  serves 
to  replace  each  pk  with  a  new  perturbed  by  some  KAM  displacement  angle  dpk-  The 
displacement  manifests  in  the  KAM  series  as 

g  (0  =  ^0  +  ^  [C-jCos  {j  ■  {p  +  5p})  +  S-jSin  (j  ■  {p  +  (5(^})] 
j 

=  ^0  +  ^  [C-jCos  {j  ■p  +  j  ■  Sp)  +  S-jSin  {j  ■  p  +  j  ■  5p)\  (3.28) 

j 

Eqn.  (3.28)  describes  in  general  terms  a  new  state  at  time  t  with  the  aforementioned 
displacement.  What  we  usually  seek  in  formation  design,  however,  is  initial  conditions 
for  two  satellites  in  physical  coordinates,  {q,p),  which  will  allow  for  the  desired  relative 
motion.  We  may  use  Eqn.  (3.6)  to  find  the  initial  coordinates  for  some  Satellite  1  at  an 
initial  time  to  =  0  simply  by: 

QSiito)  =  Co +  '^Cj  (3.29) 

3 

Given  some  initial  desired  KAM  angular  separation  vector  5p  between  the  satellites,  the 
initial  coordinates  for  Satellite  2  may  be  found  from  Eqn.  (3.28)  as 

qs2{to)  =  Co  +  ^  [C-jCos  {j  ■  5p)  +  S-jSin  (j  •  dp)]  (3.30) 

3 

In  order  to  find  the  initial  values  of  the  momenta  for  satellites  1  and  2  from  the  series,  we 
realize  that,  for  a  satellite  Sn  in  the  ECEE  frame  [25], 

PS„(C  =  T^2  PsSt)  +  qsnit)  (3.31) 

where  7^2  is  the  identity  matrix  I3X3  and 

0  — 6^0  0 

=  w©  0  0  (3.32) 

0  0  0 
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The  values  of  qs„{t)  in  Eqn.  (3.31)  may  be  found  from  the  series  (3.6)  by  simply  differen¬ 
tiating  with  respect  to  time,  so  that 

^ (t)  =  ^  [-{j  ■  n)CjSin  {j  ■  nt)  +  {j  ■  n)SjCos  {j  ■  ot)]  (3.33) 

] 

For  a  description  of  the  implementation  of  the  above  process,  continue  to  the  next  chapter. 


35 


IV.  Results  and  Findings 

This  chapter  will  present  some  of  the  results  of  the  current  work.  It  will  begin  with  the 
outcomes  of  the  KAM  torus  fitting  process  for  various  orbits,  including  an  analysis  of 
accuracy  gains  as  a  function  of  the  number  of  KAM  series  terms  included.  Next,  the 
nature  of  satellite  separations  in  the  torus  angles  ipk  will  be  discussed,  followed  by  the 
results  of  the  actual  satellite  formation  analyses.  Finally,  a  brief  sketch  on  the  effects  on 
KAM  tori  of  non-conservative  forces  will  be  given. 

4-1  KAM  Torus  Fitting 

This  section  shows  results  from  the  process  of  fitting  KAM  tori  to  integrated  data 
of  four  satellite  orbits  which  are  combinations  of  two  altitudes  (320  km  and  630  km) 
and  two  inclinations  (i  =  15°,  30°,)  .  The  orbits  were  chosen  as  representative  orbits  of 
common  Earth  satellites  and  serve  to  demonstrate  the  method  in  question.  The  method 
may,  of  course,  be  generalized  (within  certain  bounds)  to  obtain  KAM  tori  for  a  much 
greater  range  of  orbit  altitudes  and  inclinations.  The  intricacies  of  the  KAM  tori  for  the 
examined  orbits  show  definite  patterns,  thus  allowing  the  analyst  to  make  useful  statements 
about  the  general  case.  It  is,  however,  beyond  the  scope  of  this  work  to  complete  a  full 
survey  of  KAM  tori  for  all  possible  inclination/altitude  combinations,  as  the  fitting  and 
decomposition  process  can  be  quite  time-consuming  and  computationally  intensive. 

4.1.1  Torus  series  terms  and  aeeuraey  gains.  When  fitting  KAM  tori  to  orbital 
data,  the  issue  of  how  far  to  extend  the  approximation  naturally  arises.  While,  in  general, 
including  more  terms  in  the  series  given  by  Eqn.  (3.6)  tends  to  increase  the  accuracy  of 
results,  a  point  is  reached  similar  to  that  encountered  in  other  approximation  analyses 
where  the  dynamicist  encounters  the  law  of  diminishing  returns  -  that  is,  there  is  a  certain 
number  of  terms  for  each  orbit  after  which  gains  in  accuracy  come  extremely  slowly,  if  at 
all.  In  an  effort  to  numerically  characterize  the  relationship  between  the  number  of  Fourier 
series  terms  and  the  accuracy  of  an  orbit  reconstructed  from  the  series,  the  following 
procedure  was  utilized.  First,  the  orbit  for  a  satellite  at  a  specific  altitude/inclination 
combination  was  numerically  integrated  to  obtain  approximately  12  months  of  orbital  data 
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with  a  windowed  integration  error  similar  in  order  to  that  shown  in  Figure  (3.8).  Then, 
the  basis  frequencies  are  determined  to  a  very  high  accuracy  using  an  appropriate  window 
function,  as  in  Chapter  III.  Next,  a  first-pass  KAM  torus  approximation  was  calculated  by 
using  Laskar’s  and  Wiesel’s  methods  of  Fourier  analysis,  using  a  certain  index  limit  vector. 
This  approximation  gives  us  a  good  starting-point  for  the  future  spectral  decomposition,  as 
it  identifies  the  locations  of  the  highest  peaks.  Also,  one  may  use  this  torus  approximation 
as  a  check  to  determine  whether  or  not  the  index  limits  are  sufficient  to  capture  the  spectral 
behavior  of  the  dynamics:  finding  the  power  spectral  density  plot  of  the  orbit  “calculated” 
from  the  approximate  series  and  then  comparing  that  PSD  to  the  PSD  calculated  from 
the  “true” ,  integrated  orbital  data  supplies  the  analyst  with  critical  knowledge  of  missing 
peaks.  For  an  example  of  this,  see  Fig.  (4.1)  below;  in  this  figure,  the  red  curve  represents 
the  partial  PSD  of  the  integrated  ( “true” )  data  and  the  blue  curve  represents  the  partial 
PSD  of  the  orbit  constructed  using  the  fitted  torus  series  with  index  limits  jum  =  [5,5, 1]. 
The  black  diamonds  in  the  figure  show  the  locations  of  the  frequency  combinations  included 
in  the  analysis  -  i.e.,  combinations  of  the  basis  frequencies  up  to  the  limits  in  jum-  In  this 
case,  the  basis  frequency  set  was  determined  to  be  D  ss  [0.93243, 0.060025,  0.0018922].  The 
figure  clearly  shows  that  the  largest  peaks  in  this  partial  frequency  window  occur  at  the 
combinations  of  (e.g.  the  peak  at  cu  ~  0.874  is  the  combination  cj  =  Di  —  D2,  etc.).  It  is 
also  clear  that,  due  to  the  limitations  of  the  series  indices  included,  there  are  peaks  that  are 
“missed”  in  the  analysis.  Please  note  that  this  jum  is  merely  an  example  case  and  does  not 
represent  a  typical  index  limit  to  obtain  a  highly  accurate  KAM  torus.  For  informational 
purposes,  the  coordinate  residuals  for  the  approximate  torus  determined  with  these  index 
limits  as  compared  to  the  integrated  data  are  given  in  Fig.  (4.2)  below. 

After  finding  the  approximate  torus  series,  we  may  input  said  calculated  series  along 
with  the  original  orbital  data  into  the  spectral  decomposition  routine,  where  the  signal 
is  decomposed  into  its  true  frequency  contributions  up  to  some  maximum  series  term 
limit.  The  resulting  “refined”  series  provides  a  much  higher  accuracy  than  the  approximate 
Fourier  analysis  can,  and  we  may  use  the  series  to  investigate  the  title  question  of  this 
subsection:  how  the  accuracy  depends  on  the  number  of  series  terms.  We  proceed  to 
reconstruct  orbital  data  from  the  KAM  series  by  calculating  the  coordinates  and  momenta 
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Figure  4.1. 


Figure  4.2. 


Partial  PSDs  of  320  km,  i  =  30°  orbit  (red)  and  approximate  series  with 
jiim  =  [5,5, 1]  (blue) 
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at  every  time  t  using  Eqns.  (3.6),  (3.31)  and  (3.33)  and  the  Hfc,  Cj  and  Sj  from  the  refined 
series.  These  calculations  may  be  conveniently  performed  to  include  any  number  of  the 
terms  analyzed  in  the  spectral  decomposition;  the  results  of  the  12-month  RMS  residuals 
for  each  coordinate  are  shown  versus  the  number  of  included  series  terms  in  Fig.  (4.3) 
below.  The  orbit  is  the  same  320km,  30°  orbit  discussed  above,  but  the  KAM  fitting  was 


Figure  4.3.  Coordinate  RMS  residuals  over  1  year  vs.  number  of  series  terms  for  320km, 
30°  orbit 

performed  with  jum  =  [20,  25, 2]  to  obtain  frequency  combinations  up  to  the  Nyquist  limit. 
To  be  clear  regarding  convention,  the  “number  of  series  terms”  in  the  figures  below  refers 
not  to  the  jum,  but  rather  is  related  to  the  number  of  terms  included  from  the  ordered  set 
W  of  frequency  combinations  discussed  in  §3. 1.2.4.  The  “number  of  series  terms”  in  the 
plots  specifically  refers  to  the  number  of  C  and  S  coefficients  used  and  is  calculated  from 
the  number  of  frequency  combinations  as 

^terms  —  2A(^  (4:-l) 

since  there  are  two  trigonometric  coefficients  for  each  frequency  combination  (i.e.  a  real 
and  imaginary  part  of  each  complex  amplitude).  The  RMS  residuals  for  each  coordinate 
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are  determined  in  the  usual  fashion  as 


rRMSk 


'i 


1 

M 


M 

^  ^  Qk{tm) 

m=l 


(4.2) 


Similarly,  the  RMS  errors  in  the  momenta  for  the  aforementioned  320km,  30  degree  orbit 
can  be  seen  in  Figure  (4.4).  It  is  easily  seen  from  Figs.  (4.3)  and  (4.4)  that,  for  the 


Figure  4.4.  Momentum  RMS  residuals  over  1  year  vs.  number  of  series  terms  for  320km, 
30°  orbit 

orbit  examined,  there  is  little  relative  accuracy  gain  after  approximately  1500  terms  in  the 
momenta  and  after  approximately  2000  terms  in  the  coordinates.  It  is  instructive  to  note 
that,  as  mentioned  previously,  1500  terms  correspond  to  roughly  N/2  or  750  frequency 
combinations;  similarly,  2000  terms  correspond  to  1000  frequency  combinations. 

Extension  of  this  methodology  to  an  orbit  of  630  km  and  i  =  30°  gives  the  expected 
result:  a  higher  altitude  orbit  is  affected  less  by  the  geopotential  perturbations,  and  may 
thus  be  represented  to  the  maximum  possible  accuracy  by  fewer  terms  than  the  equivalent 
low-altitude  case.  Figures  (4.5)  and  (4.6)  below  show  the  RMS  residuals  for  the  coordinates 
and  momenta,  respectively,  for  the  630  km  orbit  over  a  period  of  1  year.  As  seen  in  the 
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figures,  the  accuracy  reaches  its  steady  state  after  approximately  1200  terms  (600  sets  of 
C  and  S  coefficients)  in  both  the  coordinates  and  momenta. 


Figure  4.5.  Coordinate  RMS  residuals  over  1  year  vs.  number  of  series  terms  for  630km, 
30°  orbit 

It  is  important  to  note  that,  in  the  torus  construction  process  utilized  in  this  work, 
the  accuracy  of  the  calculated  momenta  is  always  lower  than  that  of  their  associated 
coordinates.  This  is  due  to  the  fact  that  the  error  in  the  momenta  is  a  compound  error 
composed  of  the  errors  of  the  coordinates  and  the  coordinate’s  derivatives  q^,  as  is 
obvious  from  cursory  examination  of  Equations  (3.31)  and  (3.33). 
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Figure  4.6.  Momentum  RMS  residuals  over  1  year  vs.  number  of  series  terms  for  630km, 
30°  orbit 

4-1-2  Orbit  1:  Nearly  circular,  320km  Altitude  and  i  =  30°.  The  first  orbit, 
involving  an  altitude  of  320km  and  an  inclination  of  30°,  is  mentioned  extensively  in  the 
previous  subsection  and  is  used  in  this  work  as  a  sort  of  base  orbit.  As  stated  above, 
the  satellite’s  motion  for  orbit  1  was  integrated  per  the  methodology  in  Chapter  III  over 
a  period  of  time  [— T,  T]  where  T  ~  6  months,  yielding  a  total  integration  period  of 
approximately  twelve  months,  or  one  year.  Following  the  integration,  a  first-pass  torus 
was  constructed  and  then  the  data  were  spectrally  decomposed  using  the  method  out¬ 
lined  in  Chapter  III  and  §4.1.1  above.  The  residuals  were  then  found  by  subtracting 
the  “calculated”  steps  using  the  KAM  series  at  each  time  step  from  the  “true”  states  of 
the  integrated  orbit.  Figures  (4.7)  and  (4.8)  below  show  the  coordinate  and  momentum 
residuals  for  the  320km,  i  =  30°  orbit  before  and  after  decomposition;  i.e.  with  only  the 
first-pass  torus  approximation  (Fig.  (4.7))  and  with  the  more  accurate  decomposed  torus 
(Fig.  (4.8)).  The  torus  used  to  obtain  the  residuals  in  the  figures  uses  maximum  index 
limits  jiim  =  [20,25,2].  We  note  that  the  q2  residuals  in  Fig  (4.7)  are  not  visible  because 
they  are  practically  equivalent  to  and,  therefore,  occluded  by  the  qi  curve. 
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Figure  4.7.  Coordinate  and  momentum  residuals  over  1  year  before  spectral  decompo¬ 
sition  for  nearly  circular,  320km,  30°  orbit 


Figure  4.8.  Coordinate  and  momentum  residuals  over  1  year  after  spectral  decomposi¬ 
tion  for  nearly  circular,  320km,  30°  orbit 


4-. 1.3  Orbit  2:  Nearly  circular,  320km  Altitude  and  i  =  15°.  An  orbit  at  the 
same  320km  altitude  but  at  i  =  15°  was  analyzed  to  obtain  the  coordinate  and  momentum 
residuals  in  the  below  Fig.  (4.9).  With  the  lower  inclination,  the  orbit  proved  to  be 
much  more  easily  analyzed.  The  coordinate  residuals  without  any  decomposition  were 
practically  identical  to  those  shown  in  the  figure;  however,  the  decomposition  improved 
the  momentum  residuals  by  several  orders  of  magnitude.  The  series  was  decomposed  to 
only  1400  terms  in  this  case.  There  is  an  evident  demonstrated  trend  towards  generally 
easier  calculation  of  tori  with  lower  inclinations  -  see  §5.1.1  Limitations  and  considerations 
in  the  KAM  fitting  process  in  Chapter  V. 
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Figure  4.9.  Coordinate  and  momentum  residuals  over  1  year  after  spectral  decomposi¬ 
tion  for  nearly  circular,  320km,  15°  orbit 


4-1-4  Orbit  3:  Nearly  circular,  630km  Altitude  and  i  =  30°.  The  third  orbit 
used  throughout  the  rest  of  this  work  is  an  orbit  at  an  altitude  of  630km  and  i  =  30°, 
and  its  residuals  are  displayed  in  Fig.  (4.10).  The  residuals  are  generally  about  1/3  the 
magnitude  of  those  of  the  320km  orbit  at  the  same  inclination.  It  is  the  opinion  of  the 
author  that  this  is  a  result  of  the  orbit’s  higher  altitude,  which  allows  it  to  experience  less 
of  the  perturbation  from  the  geopotential. 
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Figure  4.10.  Coordinate  and  momentum  residuals  over  1  year  after  spectral  decomposi¬ 
tion  for  nearly  circular,  630km,  30°  orbit 


4-1.5  Orbit  4-  Nearly  circular,  630km  Altitude  and  i  =  15°.  The  final  orbit,  at 
630km  and  i  =  15°,  shows  an  unexpected  result  in  that  the  residuals  (shown  in  Fig.  (4.11)) 
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are  along  the  same  order  as  those  of  the  higher  inclination  orbit  at  the  same  altitude.  It  is 
believed  that  the  potential  perturbations  are  just  subtle  enough  at  630km  that  the  KAM 
torus  construction  process  is  limited  by  the  computational  accuracy  of  the  integration, 
which  would  account  for  the  similar  magnitudes  between  the  orbits  at  the  two  different 
inclinations. 


Figure  4.11.  Coordinate  and  momentum  residuals  over  1  year  after  spectral  decomposi¬ 
tion  for  nearly  circular,  630km,  15°  orbit 


4-1.6  Torus  actions  and  their  constancy.  Since  we  have,  by  following  the  method¬ 
ology  outlined  in  this  paper,  ostensibly  constructed  tori  in  the  Action-angle  space  (/,¥?), 
a  logical  check  is  to  calculate  and  examine  the  torus  actions  I.  As  mentioned  in  Chapter 
III,  the  torus  actions  should  be  constant  for  a  system  with  only  conservative  perturba¬ 
tions.  When  one  has  in  hand  the  torus  parameters,  the  actions  may  be  found  using  the 
Hamilton- Jacobi  theorem  (see  e.g.  [25],  [14])  by  the  contour  integral 

^  P  •  (4.3) 
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where  p  is  the  vector  of  the  momenta  as  found  by  Eqn.  (3.31)  and  the  differential  dq  for 
contour  Tj  is  found  by 

uQ  — 

d(pi 

=  [~0)CjSin  {j  ■  (pt)  +  {j)SjCos  {j  ■  (fit)]  (4.4) 

] 

To  check  the  constancy  of  the  actions,  the  above  integral  was  taken  over  the  full  27r  range 
of  (fi  and  (p2  for  the  320km,  30°  orbit;  specifically,  the  action  was  calculated  using  Eq. 
(4.4)  around  each  irreducible  contour  Ej  at  20  equally  spaced  locations  around  the  torus 
(in  intervals  of  vr/lO  from  (pi,(p2  =  [0,27r]).  Eigure  (4.12)  below  shows  the  result  of  the 
action  calculations  at  these  locations.  The  plot  clearly  shows  that  the  torus  actions  are 
constant  to  within  computational  uncertainty.  The  graphical  representation  also  shows  the 
proportions  of  the  actionsf:  the  two  actions  corresponding  to  the  angles  pi ,  p2  are  four  to 
five  orders  of  magnitude  greater  than  the  action  corresponding  to  Indeed,  these  are 
the  same  action  proportions  used  to  construct  the  3-toroidal  “projection”  into  2-toroidal 
space,  given  in  Eigure  (2.7)  in  §2.3.1. 
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Eigure  4.12.  Calculated  torus  actions  at  varying  locations  in  pi  and  p2  for  320km,  30° 
orbit 
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4-^  Satellite  formations  as  KAM  torus  angle  displacements 

This  section  will  give  the  nature  and  results  of  formation  flight  on  the  KAM  tori. 
First,  a  description  of  the  interpretation  of  torus  angle  displacements  will  be  given,  followed 
by  the  results  of  a  survey  of  satellite  drift  after  such  separations.  Then,  the  results  of  a 
more  focused  study  on  tight  formations  will  be  discussed. 

4-2.1  Nature  of  ipk  displacements.  To  create  a  formation  of  two  satellites,  we 
choose  a  reference  satellite  in  an  orbit  which  lies  on  the  torus  and  then  create  a  second 
satellite  whose  orbit  also  lies  on  the  torus,  but  separate  it  from  the  reference  satellite  by  a 
change  in  one  or  more  of  the  torus  angles  (ph,  per  §3.3.  Since,  as  mentioned  in  the  section 
entitled  On  the  visualization  of  tori,  the  torus  space  for  3-tori  such  as  these  are  not  easily 
intuitive,  it  is  instructive  to  investigate  the  manifestations  of  torus-surface  displacements  in 
cartesian  space.  The  following  subsections  give  a  short  exposition  of  how  these  separations 
in  the  pk  may,  when  possible,  be  interpreted  in  a  physical  sense. 

4-2. 1.1  Displacements  in  p\.  Displacements  in  the  torus  angle  pi  are 
perhaps  the  most  intuitive  and  correspond  to  separations  in  the  largest  basis  frequency, 
Di.  Following  the  discussion  in  Chapter  III  regarding  the  interpretations  of  the  basis 
frequencies  (c.f.  §3.1.3),  this  pi  separation  is,  in  physical  space,  roughly  equivalent  to  a 
separation  along  the  orbital  path  in  the  plane  of  the  orbit.  Figure  (4.13)  below  displays  an 
example  of  displacement  along  pi.  In  the  figure,  the  coordinate  frame  is  the  Earth-centered 
inertial  frame.  The  blue  X  represents  the  starting  position  of  the  reference  “chief”  satellite, 
the  red  X  is  the  starting  position  of  the  displaced  “deputy”  satellite,  and  the  blue  and  red 
lines  represent  the  trajectories  of  the  reference  and  displaced  satellites,  respectively.  For 
demonstration  purposes,  the  two  satellites  are  propagated  for  less  than  one  orbit  to  clearly 
show  the  initial  separation. 

4-2. 1.2  Displacements  in  p^.  Separations  purely  in  the  second  torus  angle, 
P2,  manifest  in  the  physical  space  as  a  rotation  of  the  orbit  plane  itself  around  the  ECI 
3-axis,  concordant  with  the  discussion  in  §3.1.3;  in  other  words,  inducing  a  separation  in  p2 
is  to  effectively  “force”  the  regression  of  the  node  through  a  certain  amount  before  placing 
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Figure  4.13.  Inertial  trajectories  of  satellites  after  separation  of  27r/3  in  (pi 

the  second  satellite.  Figure  (4.14)  shows  3D  plots  of  such  a  separation  of  magnitude  pi/ A 
in  the  ECI  frame:  the  left  plot  of  Fig  (4.14)  gives  a  perspective  view,  while  the  right  plot 
displays  a  “top-down”  view,  in  which  the  angular  separation  of  the  lines  of  nodes  for  the 
two  satellites  is  obvious.  (Note:  Although  the  orbits  appear  elliptical  in  the  top  view  in 
Fig  (4.14),  they  are  in  fact  circular  and  inclined,  which  creates  an  elliptical  projection  on 
the  ECI  X-Y  plane.)  It  is  important  to  note  that,  regardless  of  the  separation  amount,  the 
satellites’  trajectories  still  lie  on  the  same  torus  if  they  are  created  using  the  same  KAM 
series. 


Eigure  4.14.  Two  views  of  inertial  trajectories  of  satellites  after  separation  of  pi /A  in  ip2 
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4-2. 1.3  Displacements  in  (p^.  The  third  torus  angle  is  arguably  the  most 
difficult  to  intuit  and,  since  the  magnitude  of  the  third  fundamental  frequency  is  so 
small,  is  the  most  susceptible  to  modeling  errors  and  is  hardest  to  account  for.  Recall  from 
§3.1.3  that  the  third  fundamental  frequency  in  our  case,  fls,  is  representative  of  a  motion 
akin  to  the  apsidal  regression  rate.  Since  the  orbits  considered  in  this  work  are  nearly,  but 
not  exactly,  circular  (meaning  there  is  still  some  effective  eccentricity),  rotation  through 
the  third  torus  angle  have  the  effect  in  physical  space  of  displacement  in  the  orbit’s  plane 
in  a  small  ellipse  with  axes  proportional  to  the  residual  eccentricity  of  the  orbit.  Figure 
(4.15)  shows  the  inertial  positions  of  the  displacements  resulting  from  varying  over  the 
interval  [0,27r],  in  the  case  of  a  320km,  i  =  30°  orbit  where  the  eccentricity  e  ~  0.0013. 
The  blue  X  marks  represent  the  displacements  after  each  equal  increment  of  vr/lO  in 
up  to  27r.  We  note  again  that  this  figure  is  not  an  orbit  itself,  but  rather  for  displacements 
in  starting  position  from  the  reference  {ps  =  0),  marked  with  the  red  star  in  the  figure. 


Figure  4.15.  Inertial  displacements  for  =  [0,27r]  manifested  in  physical  space  for  a 
320km,  30°  orbit  with  e  ~  0.0013 


4-2.2  Formation  drift  survey  after  initial  toroidal  angular  displaeements  p^.  To 
begin  the  formation  analysis,  the  separation  drift  of  satellites  was  examined  for  varying 
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separations  purely  in  the  angles  (pk^  as  characterized  in  the  previous  section.  As  an  example 
of  the  results,  consider  the  case  case  involving  ipi  separations  in  the  320km,  30°  orbit. 
Figures  (4.16)  and  (4.17)  show  the  magnitude  of  the  satellites’  relative  distance  minus  the 
original  separation  distance  over  time  for  several  different  initial  separation  amounts.  The 
initial  separation  listed  at  the  top  of  each  subplot  is  the  direct  distance  from  one  satellite  to 
the  other  found  by  subtracting  the  position  vectors  of  the  satellites  at  time  to-this  should 
not  to  be  confused  with  the  arc  distance.  The  separation  angle  given  for  each  subplot  is 
in  terms  of  pi  angle,  which  is  what  leads  to  the  initial  separation  distance.  Again,  the 
distance  given  on  the  y-axis  in  the  figures  is  the  distance  obtained  by  taking  the  magnitude 
of  the  vector  from  one  satellite  to  the  other  at  each  time  and  then  subtracting  the  initial 
separation  vector  magnitude.  These  figures  clearly  show  a  result  which  proved  to  be 
predominant  throughout  this  work;  it  appears  that  when  satellites  are  put  into  formation 
according  to  KAM  torus  angles,  there  are  two  main  results:  an  oscillation  in  the  satellite 
separation  proportional  to  the  initial  separation  between  the  satellites  and  a  secular  drift 
to  some  extent.  In  many  cases,  especially  when  investigating  separations  in  ip2,  the  secular 
drift  is  of  such  small  magnitude  compared  to  the  magnitude  of  the  oscillation  that  it  is 
virtually  undetectable  in  the  plots. 

To  determine  the  secular  drift  in  the  satellite  separation  distance,  the  following  simple 
technique  was  used:  the  separation  data  for  each  case  (as  plotted  in  the  aforementioned 
figures)  was  fed  into  a  routine  which  calculated  the  slope  between  each  pair  of  peaks  in 
the  oscillating  data,  calculated  the  slope  between  each  pair  of  valleys,  took  the  average 
of  each  over  the  full  time  span,  and  then  averaged  the  two  slopes.  The  resulting  data 
slope  was  then  multiplied  by  the  time  span  to  determine  the  secular  drift.  As  an  example, 
for  the  320km  30°  orbit.  Figure  (4.18)  shows  the  10-day  secular  drift  as  a  function  of  the 
ipi  separation  angle.  As  mentioned  previously,  the  amount  of  secular  drift  seems  to  be 
proportional  to  the  amount  of  initial  separation.  To  gain  a  better  idea  of  the  relationship, 
the  secular  drift  as  a  percent  of  the  initial  separation  was  plotted  with  respect  to  the  initial 
<pi  angle,  shown  in  Figure  (4.19).  The  secular  drifts  for  the  pi  case  for  this  orbit  remain  less 
than  one  percent  of  the  original  separation  distance  for  the  angular  separations  examined. 
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Figure  4.16.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in 
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Figure  4.17.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (cont.) 
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Figure  4.18. 


Figure  4.19. 
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Secular  drift  between  satellites  over  10  days  vs.  initial  separations  in 
for  320km,  i  =  30°  orbit 


Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ipi  separation  for  320km,  i  =  30°  orbit 
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The  distance  results  for  each  of  the  four  orbits  and  each  are  included  in  Appendices 
A  through  D  for  reference  and,  in  the  interest  of  simplicity  and  brevity,  will  not  be  shown 
here.  However,  Tables  4.1  through  4.4  give  a  summary  of  selected  results  from  the  drift 
analysis  survey.  There  are  several  important  points  to  note  from  these  data.  First,  as 
shown  above,  the  oscillatory  amplitude  and  secular  drift  both  increase  with  increasing 
initial  separation  distance;  however,  the  drift  as  a  percentage  of  the  original  separation 
stays  relatively  low,  being  always  below  1  percent  (and  usually  below  0.5  percent)  for 
displacements  in  pi  and  p2  in  the  cases  surveyed.  Second,  as  shown  in  the  tables  and  in 
the  figures  in  the  Appendices  involving  p2^  there  is  a  seemingly  bounded  oscillation  with 
very  little  relative  secular  growth.  The  oscillation  is  unavoidable  and  is  a  direct  result  of 
the  separation  of  the  orbital  planes  by  rotation  about  the  ECI  z-axis  (recall  the  example 
separation  in  Figure  (4.14)).  Third,  drift  and  oscillation  after  separations  seem  to  be 
much  more  erratic  and  ill-behaved  than  after  separations  in  pi  and  p2-  It  is  the  speculation 
of  this  author  that,  as  alluded  to  in  §4. 2. 1.3,  the  very  small  but  nonzero  eccentricity  causes 
a  low  magnitude  oscillation  from  the  change  in  position  of  the  orbit’s  argument  of  perigee 

UJp. 
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Table  4.1.  Selected  drift  results  of  initial  displacements  in  ipk  for  320km,  15°  orbit  over 
10  days 


Torus 

angle 

Initial 

separation  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.049e+004 

27.04 

1.73 

0.01649 

1 

2.79 

3.252e+005 

838.2 

-186.2 

-0.05726 

1 

7.29 

8.493e+005 

2186 

-1007 

-0.1185 

1 

9.99 

1.163e+006 

2988 

-1243 

-0.1069 

1 

12.69 

1.476e+006 

3784 

-1100 

-0.07449 

2 

0.36 

4.209e+004 

1378 

-7.473 

-0.01775 

2 

5.76 

6.732e+005 

2.203e+004 

281.1 

0.04176 

2 

21.96 

2.552e+006 

8.358e+004 

-2869 

-0.1124 

2 

48.96 

5.552e+006 

1.818e+005 

-2833 

-0.05102 

2 

70.56 

7.738e+006 

2.534e+005 

-4991 

-0.06449 

2 

92.16 

9.65e+006 

3.161e+005 

-9286 

-0.09622 

2 

103 

1.048e+007 

3.433e+005 

-1.012e-h004 

-0.09651 

3 

0.36 

122.3 

59.52 

0.7927 

0.6483 

3 

5.76 

1941 

952 

12.56 

0.6468 

3 

21.96 

7126 

3609 

144.2 

2.023 

3 

48.96 

1.42e+004 

7851 

339.8 

2.392 

3 

70.56 

1.789e+004 

1.095e+004 

595.6 

3.329 

3 

92.16 

1.967e+004 

1.365e+004 

-962.7 

-4.895 

3 

103 

1.993e+004 

1.482e+004 

-267.3 

-1.341 
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Table  4.2.  Selected  drift  results  of  initial  displacements  in  ipk  for  320km,  30°  orbit  over 
10  days 


Torus 

angle 

Initial 

separation  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.05e+004 

19.62 

34.41 

0.3278 

1 

0.99 

1.155e+005 

215.8 

314.3 

0.2722 

1 

2.79 

3.254e+005 

608.4 

479.6 

0.1474 

1 

7.29 

8.496e+005 

1589 

-1183 

-0.1393 

1 

9.99 

1.164e+006 

2175 

-2624 

-0.2255 

1 

12.69 

1.477e+006 

2756 

-3350 

-0.2268 

2 

0.36 

4.209e+004 

5416 

221.6 

0.5265 

2 

5.76 

6.731e+005 

8.656e+004 

4463 

0.6631 

2 

21.96 

2.552e+006 

3.285e+005 

1.132e-h004 

0.4435 

2 

48.96 

5.551e+006 

7.148e+005 

2.418e-h004 

0.4356 

2 

70.56 

7.738e+006 

9.964e+005 

3.278e-h004 

0.4236 

2 

92.16 

9.65e+006 

1.243e+006 

3.977e-h004 

0.4122 

2 

103 

1.048e+007 

1.35e+006 

4.142e-h004 

0.3951 

3 

0.36 

105.6 

54.32 

-0.7253 

-0.6869 

3 

5.76 

1662 

868.6 

-39.27 

-2.363 

3 

21.96 

5942 

3292 

-14.83 

-0.2496 

3 

48.96 

1.133e+004 

7164 

-141.6 

-1.249 

3 

70.56 

1.385e+004 

9985 

209.5 

1.513 

3 

92.16 

1.502e+004 

1.246e+004 

28.76 

0.1915 

3 

103 

1.531e+004 

1.353e+004 

-657.1 

-4.291 
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Table  4.3.  Selected  drift  results  of  initial  displacements  in  ipk  for  630km,  15°  orbit  over 
10  days 


Torus 

angle 

Initial 

separation  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.098e+004 

25.34 

5.025 

0.04577 

1 

0.99 

1.208e+005 

278.8 

40.52 

0.03354 

1 

2.79 

3.404e+005 

785.5 

36.21 

0.01064 

1 

7.29 

8.888e+005 

2049 

-278.7 

-0.03135 

1 

9.99 

1.217e+006 

2801 

-468.8 

-0.03851 

1 

12.69 

1.545e+006 

3548 

-508.8 

-0.03293 

2 

0.36 

4.404e+004 

1449 

104.6 

0.2376 

2 

5.76 

7.043e+005 

2.317e+004 

2052 

0.2914 

2 

21.96 

2.67e+006 

8.788e+004 

5043 

0.1889 

2 

48.96 

5.809e+006 

1.912e+005 

1.305e-h004 

0.2246 

2 

70.56 

8.096e+006 

2.665e+005 

1.826e-h004 

0.2256 

2 

92.16 

l.Ole+007 

3.323e+005 

2.254e-h004 

0.2233 

2 

103 

1.097e+007 

3.61e+005 

2.304e-h004 

0.2101 

3 

0.36 

116.8 

57.24 

-2.388 

-2.046 

3 

5.76 

1853 

915.4 

2.153 

0.1162 

3 

21.96 

6794 

3468 

-2.579 

-0.03796 

3 

48.96 

1.351e+004 

7547 

146.1 

1.081 

3 

70.56 

1.699e+004 

1.052e+004 

553.3 

3.256 

3 

92.16 

1.866e+004 

1.312e+004 

119.3 

0.6394 

3 

103 

1.89e+004 

1.425e+004 

180.8 

0.9568 
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Table  4.4.  Selected  drift  results  of  initial  displacements  in  ipk  for  630km,  30°  orbit  over 
10  days 


Torus 

angle 

Initial 

separation  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.098e+004 

18.44 

23.72 

0.216 

1 

0.99 

1.208e+005 

202.8 

237.5 

0.1966 

1 

2.79 

3.405e+005 

571.5 

504.9 

0.1483 

1 

7.29 

8.891e+005 

1492 

60.94 

0.006854 

1 

9.99 

1.218e+006 

2043 

-651.9 

-0.05353 

1 

12.69 

1.546e+006 

2588 

-1260 

-0.08151 

2 

0.36 

4.403e+004 

5694 

-46.43 

-0.1054 

2 

5.76 

7.042e+005 

9.105e+004 

-474.3 

-0.06735 

2 

21.96 

2.67e+006 

3.456e+005 

-7336 

-0.2748 

2 

48.96 

5.808e+006 

7.517e+005 

-1.392e-h004 

-0.2397 

2 

70.56 

8.096e+006 

1.048e+006 

-1.955e-h004 

-0.2415 

2 

92.16 

l.Ole+007 

1.307e+006 

-2.451e-h004 

-0.2428 

2 

103 

1.097e+007 

1.419e+006 

-2.791e-h004 

-0.2545 

3 

0.36 

101.3 

52.6 

-1.822 

-1.798 

3 

5.76 

1594 

841.3 

32.33 

2.028 

3 

21.96 

5693 

3188 

-107.9 

-1.895 

3 

48.96 

1.084e+004 

6936 

-406.4 

-3.751 

3 

70.56 

1.323e+004 

9669 

87.26 

0.6597 

3 

92.16 

1.435e+004 

1.206e+004 

334.6 

2.332 

3 

103 

1.465e+004 

1.31e+004 

70.78 

0.4831 
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4-2.3  Tight  formation  detailed  analysis.  While  the  behavior  of  satellites  at  large 
displacements  on  the  torus  is  interesting,  the  real  promise  of  KAM  formation  design  seems 
to  lie  in  analysis  close-proximity  formation  flight.  To  further  investigate  the  utility  of 
a  KAM  approach  to  formation  design,  certain  special  cases  were  examined  using  similar 
methodology  to  that  of  the  previous  section.  First,  a  tight  formation  of  five  satellites 
flying  in  a  320km,  15°  orbit  was  analyzed.  The  satellites  were  separated  in  the  torus  space 
by  ,  6(pQ  =  0.0001°  in  the  four  cardinal  torus  1-2  plane  directions,  which  corresponds  to 
approximately  11.5  meters  of  physical  separation.  Figure  (4.20)  shows  the  separations  of 
the  satellites  both  in  the  torus  space  and  in  physical  cartesian  space.  After  defining 


Torus  Space  Physical  Space 


Figure  4.20.  Initial  position  in  torus  space  and  cartesian  space  of  satellite  cluster  for 
first  close  formation  analysis  in  320km,  15°  orbit,  6(p  =  0.0001° 

the  initial  conditions  of  the  satellite  as  such,  the  trajectory  of  each  satellite  was  obtained 
through  numerical  integration  over  a  period  of  60  days.  The  position  vector  for  each 
satellite  was  then  compared  to  that  of  the  chief  (center)  satellite  at  each  time  step  to 
determine  the  relative  motion  of  the  satellite  cluster.  The  drift  results  for  each  ancillary 
satellite  with  respect  to  the  chief  are  shown  in  Figure  (4.21),  which  shows  the  distance 
of  each  satellite  from  the  chief,  and  in  Figure  (4.22),  which  shows  the  drift  from  the 
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Figure  4.21. 


Figure  4.22. 


Cluster  distance  from  chief  satellite  for  first  close  formation  analysis  in 
320km,  15°  orbit,  dipo  =  0.0001° 
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Cluster  drift  from  initial  separations  for  first  close  formation  analysis  in 
320km,  15°  orbit,  =  0.0001° 
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original  separation  distance.  The  two  satellites  separated  purely  along  the  ip2  axis  at 
ip2  =  —0.0001°,  0.0001°  show,  as  expected,  the  largest  oscillation,  concurrent  with  the 
previous  discussion  on  the  nature  of  ip2  displacements.  In  general,  however,  the  maximum 
secular  drift  of  any  satellite  away  from  the  chief  is  approximately  0.02  meters  per  60 
days,  which  corresponds  to  an  extraordinarily  small  drift  velocity  of  approximately  3.86 
nanometers  per  second. 

Another  cluster  investigated  was  in  a  somewhat  dynamically  harsher  i  =  30°  orbit  at 
the  same  320km  altitude.  The  satellites  were  also  in  the  same  cruciform  configuration  in 
the  torus  space  as  the  previous  case,  with  5(po  =  0.0001°  as  given  on  the  left  plot  of  Figure 
(4.20),  with  the  positions  changing  in  proportion  in  the  physical  space  according  to  the 
different  orbital  configuration.  The  separation  plots  for  this  second  case  are  shown  in  Figs. 
(4.23)  and  (4.24).  The  maximum  secular  drift  experienced  by  the  satellites  in  the  cluster 
is  noticeably  higher  because  of  the  higher  inclination,  reaching  approximately  0.25m  over 
60  days,  corresponding  to  an  average  secular  drift  rate  of  approximately  48.2  nanometers 
per  second. 

Three  additional  trials  of  tight  formations  were  examined:  one  for  a  320km,  i  =  15° 
with  a  larger  separation  of  5ipQ  =  0.001°,  one  for  a  higher  630km,  i  =  15°  orbit  at  the 
small  separation  of  6(po  =  0.0001°,  and  one  for  a  630km,  i  =  30°  orbit  at  a  separation  of 
6ipo  =  0.001°.  Table  4.5  shows  a  summary  of  the  results  from  these  runs,  and  the  data 
may  be  seen  in  graphical  form  through  plots  in  Appendix  E.  The  table  shows  that  the 
highest  drift  rate  occurs  in  the  last  case  mentioned,  which  is  understandable  due  to  its 
higher  inclination  and  larger  physical  separation;  however,  the  drift  rate  is  still  only  on  the 
order  of  0.313  ^m/s. 
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Satellite  #  2 


Figure  4.23. 


Figure  4.24. 


Satellite  #  1 


Satellite  #  3 


Cluster  distance  from  chief  satellite  for  first  close  formation  analysis  in 
320km,  30°  orbit,  5ipo  =  0.0001° 


Satellite  #  1 


Satellite  #  3 


Time  (days) 


Cluster  drift  from  initial  separations  for  first  close  formation  analysis  in 
320km,  30°  orbit,  =  0.0001° 
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Table  4.5.  Results  of  tight  formation  analysis  for  various  orbits  and  separations 


Alt  =  320km,  i  =  15°,  Stpo  =  0.0001 

Sat.  No. 

Initial  sep. 
from  chief  (m) 

Avg.  Oscillation 
amplitude  (m) 

Drift  over 

60  days  (m) 

Drift 

rate  (m/s) 

1 

11.66 

0.03379 

0.01425 

2.75e-009 

2 

11.69 

0.3827 

-0.01179 

-2.275e-009 

3 

11.66 

0.03378 

0.01035 

1.996e-009 

4 

11.69 

0.3827 

-0.01087 

-2.097e-009 

Alt  =  320km,  i  =  lh° ,5ipo  =  0.001 

Sat.  No. 

Initial  sep. 
from  chief  (m) 

Avg.  Oscillation 
amplitude  (m) 

Drift  over 

60  days  (m) 

Drift 

rate  (m/s) 

1 

116.6 

0.3378 

0.1199 

2.312e-008 

2 

116.9 

3.827 

-0.1161 

-2.24e-008 

3 

116.6 

0.3378 

0.1164 

2.246e-008 

4 

116.9 

3.827 

-0.1159 

-2.235e-008 

Alt  =  320km,  i  =  30°,  Stpo  =  0.0001 

Sat.  No. 

Initial  sep. 
from  chief  (m) 

Avg.  Oscillation 
amplitude  (m) 

Drift  over 

60  days  (m) 

Drift 

rate  (m/s) 

1 

11.66 

0.03303 

0.2361 

4.555e-008 

2 

11.69 

1.499 

0.1301 

2.509e-008 

3 

11.66 

0.03303 

0.2378 

4.587e-008 

4 

11.69 

1.499 

0.128 

2.47e-008 

Alt  =  630km,  i  =  15°,  Stpo  =  0.0001 

Sat.  No. 

Initial  sep. 
from  chief  (m) 

Avg.  Oscillation 
amplitude  (m) 

Drift  over 

60  days  (m) 

Drift 

rate  (m/s) 

1 

12.2 

0.03303 

0.2361 

4.555e-008 

2 

12.23 

1.499 

0.1301 

2.509e-008 

3 

12.2 

0.03303 

0.2378 

4.587e-008 

4 

12.23 

1.499 

0.128 

2.47e-008 

Alt  =  630km,  i  =  30°,  =  0.001 

Sat.  No. 

Initial  sep. 
from  chief  (m) 

Avg.  Oscillation 
amplitude  (m) 

Drift  over 

60  days  (m) 

Drift 

rate  (m/s) 

1 

122 

0.3048 

1.626 

3.137e-007 

2 

122.3 

15.75 

1.029 

1.985e-007 

3 

122 

0.3048 

1.626 

3.136e-007 

4 

122.3 

15.75 

1.028 

1.984e-007 

63 


4.3  Atmospheric  drag  effects  on  KAM  tori 

When  one  examines  only  the  relative  motion  between  two  satellites  in  circular  orbits 
on  the  same  KAM  torus,  atmospheric  drag  effects  do  not  play  an  important  role,  as  long  as 
the  satellites  are  identical  in  drag  mass,  cross-sectional  area  and  shape  and,  if  unsymmetric, 
tend  to  fly  in  the  same  attitude-  these  are  certainly  the  case  for  many  groups  of  satellites  of 
interest  in  formation  design.  However,  when  positioning  of  the  satellites  in  an  “absolute” 
or  an  Earth-fixed  sense  is  important,  which  is  the  case  when,  for  example,  a  surveillance 
satellite  must  overfly  some  area  of  the  Earth  with  some  frequency,  atmospheric  drag  begins 
to  play  a  more  important  role  in  the  orbital  calculations.  Though  not  the  main  focus  of 
this  paper,  a  limited  study  was  performed  to  determine  the  effects  of  atmospheric  drag  on 
the  KAM  tori  of  Earth  satellites. 

Even  before  beginning  any  numerical  analysis  regarding  drag,  we  may  derive  some 
of  the  its  effects  in  general  on  the  KAM  torus  by  simple  intuition.  As  mentioned  in  §2.2.2, 
the  effect  of  drag  on  an  orbit  is  a  decrease  in  the  semi-major  axis  through  circularization 
and,  once  circularity  has  been  reached,  a  lowering  of  the  orbital  altitude  until  impact. 
Elementary  orbital  mechanics  shows  that  an  orbiting  object  whose  altitude  is  decreasing 
will  experience  an  increase  in  velocity.  We  would  expect,  then,  since  the  satellite  is  moving 
faster,  that  the  fundamental  KAM  torus  frequencies  would  increase  as  the  orbit  decays. 
However,  the  satellite’s  descent  due  to  drag  is  a  physical  manifestation  of  energy  degeneracy 
in  the  dynamical  system;  i.e.,  the  system  is  no  longer  isoenergetic.  This  means  that  the 
satellite  no  longer  lies  on  the  KAM  torus  we  have  so  carefully  characterized  and  calculated, 
since  one  of  the  fundamental  assumptions  made  in  Chapter  HI  involved  an  isoenergetic 
system  with  a  small,  conservative  perturbation  to  a  Hamiltonian  system. 

On  the  surface,  it  appears  that  atmospheric  drag  renders  KAM  representation  of 
orbits  useless.  On  the  contrary,  if  the  drag  force  can  be  characterized  well  enough,  it  should 
still  be  possible  to  gain  information  from  the  KAM  torus.  To  very  briefly  investigate  this, 
an  analysis  was  performed  in  which  a  satellite  began  in  a  320km,  30°  inclination  orbit. 
To  include  the  effects  of  air  drag,  this  author  proposes  the  use  of  a  Rayleigh  dissipative 
function  [14]  in  terms  of  the  qk  inserted  into  the  Hamiltonian,  so  that  Hamilton’s  functions 
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become 


qk  = 


5H_ 

Spk 

5H 


6Fr 

Pk  =  -T - 

OQk  oqk 


(4.5) 


where  x  is  a  coefficient  dependent  on  spacecraft  characteristics  and  the  atmospheric  density 
given  as 


X  =  P 


CdA 

m 


P 


(4.6) 


and  Fr  is  a  Rayleigh  function 


Fr 


1 

6 


{qi  +  q2  +  qs) 


3 

2 


(4.7) 


When  integrating  Eqns.  (4.5)  to  provide  trajectory  data,  the  instantaneous  atmospheric 
density  p  in  (4.6)  was  calculated  at  every  time  step  using  a  basic  exponential  atmospheric 
model  for  the  Earth  (c.f.  [19]  for  a  more  detailed  explanation).  A  more  accurate  density 
model,  if  available,  could  obviously  be  utilized  in  its  place. 

A  reference  trajectory  was  first  integrated  over  t  =  [—T,  T]  without  drag  and  a  KAM 
torus  constructed  for  it.  Next,  the  same  initial  conditions  were  integrated  forwards  with 
the  drag  acceleration  included.  Then,  the  state  vectors  of  this  “drag-perturbed”  trajectory 
at  each  time  from  1  to  64  time  units  were  extracted.  These  state  vectors  were  then  used  as 
initial  conditions  for  new  forward  and  backward  integrations,  again  over  t  =  [— T,  T],  and 
KAM  tori  constructed  for  each  trajectory.  In  this  way,  we  have  found,  in  essence,  osculating 
KAM  tori,  in  that  the  tori  are  what  would  be  obtained  if  one  could  instantaneously  “switch 
off”  the  drag  force  at  any  point  in  a  trajectory  and  analyze  the  motion  of  the  new  orbit. 

The  tori  basis  frequencies  are  necessarily  calculated  in  order  to  construct  new  KAM 
tori  for  each  state  of  the  drag-perturbed  trajectory.  The  set  of  the  basis  frequencies  for  the 
constructed  tori  then  will  show  the  trend  in  the  frequencies  as  drag  acts  upon  the  satellite. 
Eigures  (4.25)  through  (4.26)  show  the  change  in  each  of  the  three  basis  frequencies  as  found 
through  the  above  analysis  on  a  320km,  30°  orbit.  Since  this  was  purely  a  demonstrative 
investigation,  for  the  sake  of  computational  efficiency,  the  geopotential  expansion  was  only 
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included  to  order  and  degree  m,  n  =  4  rather  than  the  usual  m,  n  =  20  limit  used  elsewhere 
throughout  this  work.  It  is  clear  from  the  figures  that  the  fundamental  frequencies  display 


Figure  4.25.  Change  in  Oi  due  to  drag  in  a  320km,  30°  orbit,  including  geopotential 
expansion  to  m,  n  =  4 

the  expected  increase  due  to  the  lowering  of  the  orbit;  also,  all  three  display  the  same 
roughly  linear  behavior  for  this  relatively  short  time,  though  they  grow  at  different  rates 
proportional  to  the  initial  frequency  magnitude.  The  above  results  serve  to  illustrate  that, 
even  though  the  orbit  of  a  satellite  acted  on  by  an  energy-dissipative  force  like  drag  will 
not  stay  on  the  same  KAM  torus,  the  torus  evolves  in  a  way  consistent  with  the  expected 
and  observed  orbital  deterioriation. 
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0.060027 


Figure  4.26. 


Figure  4.27. 


time  (min) 

Change  in  ^2  due  to  drag  in  a  320km,  30°  orbit,  including  geopotential 
expansion  to  m,  n  =  4 
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Change  in  Q3  due  to  drag  in  a  320km,  30°  orbit,  including  geopotential 
expansion  to  m,  n  =  4 
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V.  Conclusions 


This  chapter  will  posit  conclusions  derived  from  the  results  of  the  work  presented  through¬ 
out  Chapter  IV.  Torus  construction,  its  limitations  and  applicability  will  be  discussed, 
followed  by  the  key  conclusions  of  formation  design  on  the  torus.  Finally,  recommenda¬ 
tions  for  future  work  will  be  given. 

5. 1  Torus  Construction 

As  seen  in  Chapter  IV,  the  KAM  tori  derived/calculated  for  this  work  have  been 
utilized  to  accurately  describe  the  positions  of  virtual  satellites  in  their  orbits  to  within  0.1 
meter  over  one  year  in  some  cases-  this  obviously  demonstrates  that  KAM  is  a  powerful 
tool  with  some  degree  of  utility  in  orbital  mechanics  and  design.  This  section  will  describe 
the  limitations  of  torus  construction  and  thereby  postulations  regarding  its  applicability. 

5.1.1  Limitations  and  considerations  in  the  KAM  fitting  process.  As  mentioned 
above,  the  methods  used  in  this  work  to  construct  KAM  tori  have  yielded  promising  and 
potentially  valuable  results;  these  must  be  tempered,  however,  by  the  reality  of  limitations 
inherent  in  the  process  and  their  associated  impacts  on  the  potential  usefulness  of  the 
system.  For  example,  the  results  of  Chapter  IV  showed  a  trend  in  which  tori  of  orbits 
of  higher  altitudes  were  determinable  to  higher  accuracy  using  fewer  KAM  series  terms 
than  their  lower-altitude  counterparts  sharing  the  same  inclination.  However,  the  opposite 
trend  was  found  with  regards  to  changing  inclination;  as  inclination  increases,  the  orbit’s 
torus  is  generally  more  difficult  to  obtain  to  any  high  degree  of  accuracy.  Indeed,  orbits  at 
both  320km  and  630km  altitudes  were  examined  by  the  current  author  whose  inclinations 
were  as  high  as  37r/4  ~  43°,  the  tori  of  which  were  only  able  to  be  calculated  with  residuals 
on  the  order  of  650km  over  one  year. 

This  author  believes  that  the  phenomenon  of  KAM  tori  becoming  more  difficult  to 
determine  accurately  with  increasing  inclination  is  directly  an  effect  of  the  shrinking  of 
the  third  basis  frequency  as  seen  in  Figure  (3.6)  as  i  ^  i* ,  where  i*  ~  63.4°  is  the  critical 
inclination  discussed  in  §3.1.3.  As  H3  falls  closer  and  closer  to  a  zero  value,  the  trajectory 
knowledge  must  be  more  and  more  accurate  over  a  longer  and  longer  time  to  accurately 
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determine  the  basis  frequencies  with  the  methodology  in  Chapter  III.  In  the  limit  where 
i  =  i*,  the  basis  frequency  set  O  =  [Hi,  H2,  H3]  would  be  said  to  be  commensurable  (i.e.,  two 
or  more  elements  of  the  set  have  a  common  divisor),  violating  the  diophantine  condition 
mentioned  in  [9]  and  leading  to  an  exacerbated  problem  of  small  divisors.  In  this  case, 
the  torus  would  be  practically  incalculable.  In  this  case,  however,  where  \i  —  i*\  ~  0, 
it  is  expected  that  the  torus  could  be  described  with  reasonable  accuracy  in  terms  of 
the  set  of  only  two  fundamental  frequencies,  H  =  [Hi,H2].  Such  an  investigation  was 
beyond  the  scope  of  this  work,  but  it  was  determined  that  the  transition  region  where 
0.0017  rad/ s  >  H3  >  0  leads  to  a  severe  decline  in  the  accuracy  of  the  results. 

Similarly,  one  also  encounters  a  commensurability  issue  when  two  of  the  are  very 
close  to  low  integer  multiples  of  each  other,  even  if  they  are  nowhere  close  to  zero-valued. 
This  would  be  most  frequently  encountered  as  a  relationship  between  the  two  largest 
frequencies,  and  would  only  be  soluble  with  trajectory  knowledge  of  greater  accuracy  to 
resolve  the  two  frequencies 

5.1.2  Applicability  of  KAM  theorem.  Of  pivotal  importance  in  the  KAM  pro¬ 
cess,  and  intimately  related  to  the  commensurability  considerations  discussed  above,  is 
the  trajectory  knowledge  of  the  orbit  whose  KAM  torus  is  sought.  In  the  current  work, 
the  trajectory  knowledge  was  obtained  through  numerical  integration,  which  is  much  more 
easily  performed  than  physical  on-orbit  determination  to  a  high  degree  of  accuracy.  In  the 
author’s  estimation,  the  main  limitation  on  the  applicability  of  KAM  torus  construction 
to-date  is  the  extremely  long  periods  of  accurate  trajectory  knowledge  required  to  create 
usable  KAM  tori.  It  is  indeed  quite  rare,  due  to  the  inherent  error  in  any  real-world  track¬ 
ing  or  position  determination  system  and  the  necessity  of  spacecraft  maneuvering,  to  know 
the  true  trajectory  of  a  satellite  well  enough  to  very  accurately  implement  KAM  theorem. 
Little  has  shown  the  results  of  such  investigations  of  actual  satellites’  tori.  Even  in  the 
“ideal”  case  of  numerically  integrated  data,  there  is  a  finite  limit  of  accuracy  (as  discussed 
in  Chapter  III  and  above)  which  casts  accuracy  bounds  on  KAM  torus  determination. 
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5.2  KAM  Formation  Flight 

The  results  of  Chapter  IV  show  that  an  accurate  KAM  torus  for  a  given  orbit  can 
yield  useful  fundamental  knowledge  of  the  orbit’s  behavior,  including  its  major  modes  of 
motion  and  their  frequencies  and  proportions.  Additionally,  it  was  shown  that  formations 
may  be  configured  on  the  KAM  torus  in  both  tight  (close-lying)  and  wide  formations  to 
obtain  physical  separations  with  varying  levels  of  relative  oscillation  and  secular  drift.  In 
the  survey  of  wide  formations,  the  satellites  generally  experienced  secular  drifts  of  less  than 
0.5  percent  of  the  original  separations  in  torus  angles  (fi  and  ip2  over  10  days,  even  with 
separations  of  as  much  as  1-5  million  meters.  It  seems  that  the  amount  of  secular  drift 
between  satellites  and  the  oscillation  distance  amplitude  is  proportional  to  the  amount  of 
separation,  which  casts  doubt  on  the  utility  of  this  method  for  designing  constellations  of 
satellites  separated  by  large  distances. 

The  KAM  method  seems  to  show  promise,  however,  in  genesis  of  satellite  clusters 
with  small  5(po.  Orbits  of  satellite  clusters  separated  by  ten  to  one  hundred  meters  ac¬ 
cording  to  KAM  tori  showed  drift  rates  in  the  nanometer  to  micrometer  range  over  60 
day  integration  windows.  If,  in  fact,  satellites  can  be  placed  on  a  KAM  torus  with  such 
angular  separations,  the  low  and  accurate  thrust  of  electric  propulsive  devices  could  easily 
overcome  such  a  drift  rate. 

5.3  Recommendations  for  Future  Study 

It  is  the  author’s  opinion  that,  in  order  to  transform  KAM  theory  into  a  viable  tool 
for  astrodynamicists,  certain  additional  areas  in  this  promising  field  should  be  investigated 
and  characterized.  First,  a  purely  informational  survey  of  torus  construction  for  different 
altitude  and  inclination  combinations  should  be  performed,  determining  the  steps  and  basis 
frequency  number /characteristics  required  for  each  type  of  orbit.  The  envisioned  result  of 
such  an  effort  would  be  a  sort  of  database  of  KAM  torus  parameters  for  Earth  orbits. 
Second,  investigation  of  the  possibility  of  applying  KAM  torus  construction  methodology 
to  orbits  of  varying  eccentricities  should  be  executed.  Third,  additional  work  should  be 
completed  with  a  focus  on  developing  new  and  more  efficient  techniques  for  determining 
torus  frequencies  and  torus  coefficients  -  in  essence,  to  search  for  a  more  direct  map 
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M  :  (q,p)  (/,  ip),  the  ideal  map  mentioned  at  the  beginning  of  Chapter  III.  Such  effort 

would  seek  to  reduce  or  eliminate  the  need  for  extensive  numerical  integration  and/or 
highly  accurate  trajectory  knowledge  over  exceedingly  long  times  in  constructing  KAM 
torus  models  for  an  orbit.  Fourth,  the  problem  of  maneuvering  onto  KAM  tori  should 
be  investigated.  The  author  is  aware  of  current  work  with  some  success  in  this  area, 
and  feels  strongly  that  this  is  a  vital  area  of  continued  research,  with  the  ideal  result  of 
allowing  for  flight  on  specified  KAM  tori  in  spite  of  the  intrinsically  limited  precision  of 
launch  and  orbital  insertion  capabilities.  Finally,  in  the  area  of  KAM  formation  design, 
it  would  be  beneficial  to  fully  and  extensively  characterize  the  correlation  between  torus 
angle  displacements  and  their  associated  displacements  in  the  physical  space,  with  a  focus 
on  the  ability  to  “fix”  satellites  into  the  desired  relative  physical  positions  while  harnessing 
the  drift  advantages  of  their  being  on  the  same  KAM  torus. 
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Appendix  A.  Formation  drift  survey  results,  320km,  i  =  15° 


Figure  Al. 


Secular  drift  between  satellites  over  10  days  vs.  initial  separation  for 
320km,  i  =  15°  orbit 


Table  Al. 


Drift  results  of  initial  displacements  in  for  320km,  15°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.049e+004 

27.04 

1.73 

0.01649 

1 

0.99 

1.154e+005 

297.4 

-11.37 

-0.009852 

1 

1.89 

2.203e+005 

567.8 

-74.86 

-0.03398 

1 

2.79 

3.252e+005 

838.2 

-186.2 

-0.05726 

1 

3.69 

4.301e+005 

1108 

-334.2 
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Figure  A2.  Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ipi  separation  for  320km,  i  =  15°  orbit 

Table  A2.  Drift  results  of  initial  displacements  in  ip2  for  320km,  15°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

2 

0.36 

4.209e+004 

1378 

-7.473 

-0.01775 

2 

5.76 

6.732e+005 

2.203e+004 

281.1 

0.04176 

2 

11.16 

1.303e+006 

4.263e+004 

699.2 

0.05367 

2 

16.56 

1.93e+006 

6.318e+004 

-1301 

-0.06745 

2 

21.96 

2.552e+006 

8.358e+004 

-2869 

-0.1124 

2 

27.36 

3.169e+006 

1.037e+005 

-1465 

-0.04625 

2 

32.76 

3.778e+006 

1.237e+005 

-1730 

-0.04579 

2 

38.16 

4.38e+006 

1.434e+005 

-5276 

-0.1205 

2 

43.56 

4.971e+006 

1.628e+005 

-6069 

-0.1221 

2 

48.96 

5.552e+006 

1.818e+005 

-2833 

-0.05102 

2 

54.36 

6.12e+006 

2.004e+005 

-2946 

-0.04814 

2 

59.76 

6.675e+006 

2.186e+005 

-6848 

-0.1026 

2 

65.16 

7.214e+006 

2.363e+005 

-7495 

-0.1039 

2 

70.56 

7.738e+006 

2.534e+005 

-4991 

-0.06449 

2 

75.96 

8.245e+006 

2.7e+005 

-5605 

-0.06799 

2 

81.36 

8.733e+006 

2.86e+005 

-9124 

-0.1045 

2 

86.76 

9.202e+006 

3.014e+005 

-1.025e-h004 

-0.1114 

2 

92.16 

9.65e+006 

3.161e+005 

-9286 

-0.09622 

2 

97.56 

1.008e+007 

3.3e+005 

-9012 

-0.08942 

2 

103 

1.048e+007 

3.433e+005 

-1.012e-h004 

-0.09651 
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Table  A3.  Drift  results  of  initial  displacements  in  for  320km,  15°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

3 

0.36 

122.3 

59.52 

0.7927 

0.6483 

3 

5.76 

1941 

952 

12.56 

0.6468 

3 

11.16 

3723 

1843 

-65.99 

-1.773 

3 

16.56 

5455 

2730 

113.1 

2.073 

3 

21.96 

7126 

3609 

144.2 

2.023 

3 

27.36 

8725 

4481 

-44.89 

-0.5146 

3 

32.76 

1.024e+004 

5344 

-390.1 

-3.81 

3 

38.16 

1.166e+004 

6196 

-715.9 

-6.138 

3 

43.56 

1.299e+004 

7034 

150.5 

1.159 

3 

48.96 

1.42e+004 

7851 

339.8 

2.392 

3 

54.36 

1.53e+004 

8653 

476.4 

3.113 

3 

59.76 

1.629e+004 

9439 

267.1 

1.64 

3 

65.16 

1.715e+004 

1.021e+004 

-226.4 

-1.32 

3 

70.56 

1.789e+004 

1.095e+004 

595.6 

3.329 

3 

75.96 

1.851e+004 

1.166e+004 

-186.2 

-1.006 

3 

81.36 

1.901e+004 

1.235e+004 

-250.6 

-1.319 

3 

86.76 

1.939e+004 

1.302e+004 

-286.7 

-1.478 

3 

92.16 

1.967e+004 

1.365e+004 

-962.7 

-4.895 

3 

97.56 

1.984e+004 

1.425e+004 

87.19 

0.4394 

3 

103 

1.993e+004 

1.482e+004 

-267.3 

-1.341 
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Figure  A3. 


Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  ipi  for  320km,  i  =  15°  orbit 
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Figure  A4.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  for  320km,  i  =  15°  orbit  (contd.) 
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days  vs.  initial  ip2  separation  for  320km,  i  =  15°  orbit 
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Figure  A7.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (^2  for  320km,  i  =  15°  orbit 
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Figure  A8.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (p2  for  320km,  i  =  15°  orbit  (contd.) 


79 


600 


Figure  A9. 


Figure  AlO. 


400 

200 


-1 000 ' - ^ ^ ^ ^ 

0  20  40  60  80  100 

Phi  separation  (deg) 


120 


Secular  drift  between  satellites  over  10  days  vs.  initial  separation  for 
320km,  i  =  15°  orbit 


Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  separation  for  320km,  i  =  15°  orbit 


80 


Deviation  (km) 


I  ° 

C  -0.02 
q 

."S  -0.04 

S  -0.06 


in 

phi^  Init 

1  0 

--  -0.2 

o  -0.4 

ro  -0.6 

■>  -0.8 

Init.  sep.  =  1 .941  km  =  5.76  deg  in  phi^ 


0  5  10 

Time  (days) 

Init.  sep.  =  3.723  km  =  1 1 .2  deg  in  phi^ 
0 


E 


0  5  10 

Time  (days) 

Init.  sep.  =  7.126  km  =  22  deg  in  phi^ 
0 


o 

ro 

■> 

Q 


0  5  10 

Time  (days) 

Init.  sep.  =  5.455  km  =  16.6  deg  in  phi^ 

Oi 
-1 
-21 

0  5  10 

Time  (days) 

Init.  sep.  =  8.725  km  =  27.4  deg  in  phi^ 


c 

-1 

q 

"cS 

-2 

■> 

CD 

-3 

Q 

Time  (days) 


Time  (days) 


Init.  sep.  =  1 0.24  km  =  32.8  deg  in  phi^  Init.  sep.  =  1 1 .66  km  =  38.2  deg  in  phi^ 


E 


c 

.O  _2 

4—*  *— 

q 

S  -4 

Q 


0  5  10 

Time  (days) 

Init.  sep.  =  12.99  km  =  43.6  deg  in  phi^ 


0  5  10 

Time  (days) 

Init.  sep.  =  14.2  km  =  49  deg  in  phi^ 


IKlUiUiiililUiUllkkllllllliaUlUililUiiiUiiiiJiXLlUUUiiiiili 


10 


lUUaiiiliii'iiUiiiitiJiLiaaailllliUililikli 


10 


Time  (days) 


Time  (days) 


Figure  All.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  for  320km,  i  =  15°  orbit 
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Figure  A12.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (/?3  for  320km,  i  =  15°  orbit  (contd.) 
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Appendix  B.  Formation  drift  survey  results,  320km,  i  =  30 


Figure  Bl. 


Secular  drift  between  satellites  over  10  days  vs.  initial  pi  separation  for 
320km,  i  =  30°  orbit 


Table  Bl. 


Drift  results  of  initial  displacements  in  pi  for  320km,  30°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.05e+004 

19.62 

34.41 

0.3278 

1 

0.99 

1.155e+005 

215.8 

314.3 

0.2722 

1 

1.89 

2.204e+005 

412.1 

466.4 

0.2116 

1 

2.79 

3.254e+005 

608.4 

479.6 

0.1474 

1 

3.69 

4.303e+005 

804.7 

354.9 

0.08249 

1 

4.59 

5.352e+005 

1001 

106 

0.01981 

1 

5.49 

6.4e+005 

1197 

-249.5 

-0.03898 

1 

6.39 

7.449e+005 

1393 

-694 

-0.09318 

1 

7.29 

8.496e+005 

1589 

-1183 

-0.1393 

1 

8.19 

9.544e+005 

1785 

-1692 

-0.1773 

1 

9.09 

1.059e+006 

1980 

-2178 

-0.2056 

1 

9.99 

1.164e+006 

2175 

-2624 

-0.2255 

1 

10.89 

1.268e+006 

2369 

-2975 

-0.2346 

1 

11.79 

1.373e+006 

2563 

-3221 

-0.2347 

1 

12.69 

1.477e+006 

2756 

-3350 

-0.2268 
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Figure  B2.  Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ipi  separation  for  320km,  i  =  30°  orbit 

Table  B2.  Drift  results  of  initial  displacements  in  ip2  for  320km,  30°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

2 

0.36 

4.209e+004 

5416 

221.6 

0.5265 

2 

5.76 

6.731e+005 

8.656e+004 

4463 

0.6631 

2 

11.16 

1.303e+006 

1.675e+005 

8481 

0.6511 

2 

16.56 

1.929e+006 

2.483e+005 

1.004e+004 

0.5206 

2 

21.96 

2.552e+006 

3.285e+005 

1.132e+004 

0.4435 

2 

27.36 

3.168e+006 

4.079e+005 

1.391e+004 

0.439 

2 

32.76 

3.778e+006 

4.865e+005 

1.532e+004 

0.4055 

2 

38.16 

4.379e+006 

5.641e+005 

1.593e+004 

0.3637 

2 

43.56 

4.971e+006 

6.402e+005 

1.936e+004 

0.3896 

2 

48.96 

5.551e+006 

7.148e+005 

2.418e+004 

0.4356 

2 

54.36 

6.12e+006 

7.88e+005 

2.66e+004 

0.4347 

2 

59.76 

6.674e+006 

8.594e+005 

2.829e+004 

0.4238 

2 

65.16 

7.214e+006 

9.289e+005 

3.129e+004 

0.4338 

2 

70.56 

7.738e+006 

9.964e+005 

3.278e+004 

0.4236 

2 

75.96 

8.244e+006 

1.062e+006 

3.179e+004 

0.3856 

2 

81.36 

8.732e+006 

1.125e+006 

3.234e+004 

0.3703 

2 

86.76 

9.201e+006 

1.185e+006 

3.606e+004 

0.3919 

2 

92.16 

9.65e+006 

1.243e+006 

3.977e+004 

0.4122 

2 

97.56 

1.008e+007 

1.298e+006 

4.122e+004 

0.4091 

2 

103 

1.048e+007 

1.35e+006 

4.142e+004 

0.3951 
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Table  B3.  Drift  results  of  initial  displacements  in  for  320km,  30°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

3 

0.36 

105.6 

54.32 

-0.7253 

-0.6869 

3 

5.76 

1662 

868.6 

-39.27 

-2.363 

3 

11.16 

3159 

1681 

-151.9 

-4.807 

3 

16.56 

4589 

2490 

40.08 

0.8733 

3 

21.96 

5942 

3292 

-14.83 

-0.2496 

3 

27.36 

7211 

4088 

253.3 

3.513 

3 

32.76 

8389 

4875 

163.7 

1.952 

3 

38.16 

9471 

5652 

-89.62 

-0.9463 

3 

43.56 

1.045e+004 

6415 

298.1 

2.852 

3 

48.96 

1.133e+004 

7164 

-141.6 

-1.249 

3 

54.36 

1.211e+004 

7898 

-51.77 

-0.4274 

3 

59.76 

1.279e+004 

8613 

-410.4 

-3.209 

3 

65.16 

1.337e+004 

9311 

-566 

-4.235 

3 

70.56 

1.385e+004 

9985 

209.5 

1.513 

3 

75.96 

1.425e+004 

1.064e+004 

36.09 

0.2533 

3 

81.36 

1.457e+004 

1.127e+004 

850.5 

5.839 

3 

86.76 

1.482e+004 

1.187e+004 

316.6 

2.137 

3 

92.16 

1.502e+004 

1.246e+004 

28.76 

0.1915 

3 

97.56 

1.518e+004 

1.301e+004 

347.9 

2.293 

3 

103 

1.531e+004 

1.353e+004 

-657.1 

-4.291 
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Figure  B3.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  for  320km,  i  =  30°  orbit 
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Figure  B4.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (pi  for  320km,  i  =  30°  orbit  (contd.) 
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Figure  B5. 


Figure  B6. 


Secular  drift  between  satellites  over  10  days  vs.  initial  ip2  separation  for 
320km,  i  =  30°  orbit 


Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ip2  separation  for  320km,  i  =  30°  orbit 
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Figure  B7.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (p2  for  320km,  i  =  30°  orbit 
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Figure  B8.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?2  for  320km,  i  =  30°  orbit  (contd.) 
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Figure  Bll.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (^3  for  320km,  i  =  30°  orbit 
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Figure  B12.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?3  for  320km,  i  =  30°  orbit  (contd.) 
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Appendix  C.  Formation  drift  survey  results,  630km,  i  =  15° 


Figure  Cl.  Secular  drift  between  satellites  over  10  days  vs.  initial  pi  separation  for 
630km,  i  =  15°  orbit 


Table  Cl.  Drift  results  of  initial  displacements  in  pi  for  630km,  15°  orbit  over  10  days 
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Initial 
sep.  (deg) 
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separation  (m) 
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amplitude  (m) 
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Figure  C2.  Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ipi  separation  for  630km,  i  =  15°  orbit 

Table  C2.  Drift  results  of  initial  displacements  in  ip2  for  630km,  15°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 
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0.36 

4.404e+004 
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Table  C3.  Drift  results  of  initial  displacements  in  for  630km,  15°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 
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0.36 
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57.24 
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Figure  C3.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (pi  for  630km,  i  =  15°  orbit 
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Figure  C4.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  for  630km,  i  =  15°  orbit  (contd.) 
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Figure  C5. 


Figure  C6. 


Secular  drift  between  satellites  over  10  days  vs.  initial  ip2  separation  for 
630km,  i  =  15°  orbit 


Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ip2  separation  for  630km,  i  =  15°  orbit 
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Figure  C7.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  ip2  for  630km,  i  =  15°  orbit 
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Figure  C8.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (p2  for  320km,  630km,  i  =  15°  (contd.) 
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Figure  Cll.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?3  for  630km,  i  =  15°  orbit 
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Figure  C12.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?3  for  630km,  i  =  15°  orbit  (contd.) 
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Appendix  D.  Formation  drift  survey  results,  630km,  i  =  30' 


Figure  Dl. 


Secular  drift  between  satellites  over  10  days  vs.  initial  ipi  separation  for 
630km,  i  =  30°  orbit 


Table  Dl. 


Drift  results  of  initial  displacements  in  pi  for  630km,  30°  orbit 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

1 

0.09 

1.098e+004 

18.44 

23.72 

0.216 

1 

0.99 

1.208e+005 

202.8 

237.5 

0.1966 

1 

1.89 

2.307e+005 

387.2 

401 

0.1739 
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2.79 

3.405e+005 

571.5 

504.9 

0.1483 
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4.503e+005 

755.8 

542.9 

0.1206 
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4.59 

5.601e+005 

940.1 
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0.09203 
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1124 
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0.0625 
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7.795e+005 

1308 

265.3 

0.03403 
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7.29 

8.891e+005 

1492 

60.94 

0.006854 
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8.19 

9.987e+005 

1676 

-164 

-0.01643 

1 

9.09 

1.108e+006 

1860 

-410.1 

-0.03701 
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1.218e+006 

2043 

-651.9 

-0.05353 
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10.89 
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11.79 
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Figure  D2.  Secular  drift  (in  percentage  of  initial  separation)  between  satellites  over  10 
days  vs.  initial  ipi  separation  for  630km,  i  =  30°  orbit 

Table  D2.  Drift  results  of  initial  displacements  in  <^2  for  630km,  30°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

2 

0.36 

4.403e+004 

5694 

-46.43 

-0.1054 

2 

5.76 

7.042e+005 

9.105e+004 

-474.3 

-0.06735 

2 

11.16 

1.363e+006 

1.763e+005 

-1955 

-0.1434 

2 

16.56 

2.019e+006 

2.612e+005 

-5070 

-0.2512 

2 

21.96 

2.67e+006 

3.456e+005 

-7336 

-0.2748 

2 

27.36 

3.315e+006 

4.29e+005 

-7749 

-0.2338 

2 

32.76 

3.953e+006 

5.115e+005 

-8775 

-0.222 

2 

38.16 

4.582e+006 

5.93e+005 

-1.133e-h004 

-0.2473 

2 

43.56 

5.201e+006 

6.731e+005 

-1.326e-h004 

-0.2549 

2 

48.96 

5.808e+006 

7.517e+005 

-1.392e-h004 

-0.2397 

2 

54.36 

6.403e+006 

8.286e+005 

-1.479e-h004 

-0.2311 
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59.76 

6.983e+006 

9.037e+005 

-1.615e-h004 

-0.2313 

2 

65.16 

7.548e+006 

9.768e+005 

-1.758e-h004 

-0.2329 
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70.56 

8.096e+006 

1.048e+006 

-1.955e-h004 

-0.2415 

2 

75.96 

8.626e+006 

1.116e+006 

-2.179e-h004 

-0.2526 

2 

81.36 

9.136e+006 

1.182e+006 

-2.301e-h004 

-0.2519 

2 

86.76 

9.627e+006 

1.246e+006 

-2.346e-h004 

-0.2437 

2 

92.16 

l.Ole+007 

1.307e+006 

-2.451e-h004 

-0.2428 

2 

97.56 

1.054e+007 

1.364e+006 

-2.633e-h004 

-0.2497 

2 
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Table  D3.  Drift  results  of  initial  displacements  in  (^3  for  630km,  30°  orbit  over  10  days 


Torus 

angle 

Initial 
sep.  (deg) 

Initial 

separation  (m) 

Oscillation 
amplitude  (m) 

Drift  (m) 

Drift 

(percent) 

3 

0.36 

101.3 

52.6 

-1.822 

-1.798 

3 

5.76 

1594 

841.3 

32.33 

2.028 

3 

11.16 

3030 

1628 

84.19 

2.779 

3 

16.56 

4399 

2411 

64.15 

1.459 

3 

21.96 

5693 

3188 

-107.9 

-1.895 

3 

27.36 

6905 

3960 

-116 

-1.679 

3 

32.76 

8030 

4722 

192.3 

2.395 

3 

38.16 

9062 

5472 

284 

3.134 

3 

43.56 

9998 

6210 

290.5 

2.906 

3 

48.96 

1.084e+004 

6936 

-406.4 

-3.751 

3 

54.36 

1.158e+004 

7648 

-384.6 

-3.323 

3 

59.76 

1.222e+004 

8342 

-3.641 

-0.0298 

3 

65.16 

1.277e+004 

9016 

209.4 

1.64 

3 

70.56 

1.323e+004 

9669 

87.26 

0.6597 

3 

75.96 

1.36e+004 

1.03e+004 

297.4 

2.186 

3 

81.36 

1.391e+004 

1.091e+004 

-90.39 

-0.6499 

3 

86.76 

1.415e+004 

1.15e+004 

64.09 

0.4529 

3 

92.16 

1.435e+004 

1.206e+004 

334.6 

2.332 

3 

97.56 

1.451e+004 

1.259e+004 

-271.5 

-1.872 

3 
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Figure  D3.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (/?i  for  630km,  i  =  30°  orbit 
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Figure  D4.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (pi  for  630km,  i  =  30°  orbit  (contd.) 
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Figure  D5. 


Figure  D6. 
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Secular  drift  between  satellites  over  10  days  vs.  initial  ip2  separation  for 
630km,  i  =  30°  orbit 
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days  vs.  initial  ip2  separation  for  630km,  i  =  15°  orbit 
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Figure  D7.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  ip2  for  630km,  i  =  30°  orbit 
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Figure  D8.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  (p2  for  320km,  630km,  i  =  30°  (contd.) 
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Figure  Dll.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?3  for  630km,  i  =  30°  orbit 
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Figure  D12.  Satellite  separation  deviation  from  initial  value  over  10  days  after  varying 
initial  separations  in  y?3  for  630km,  i  =  30°  orbit  (contd.) 
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Appendix  E.  Tight  formation  analysis  results  for  various  orbits 
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Figure  E2. 


Initial  position  in  torus  space  and  inertial  cartesian  space  of  satellite  cluster 
for  tight  formation  analysis  in  320km,  15°  orbit,  dp  =  0.0001° 


Cluster  distance  from  chief  satellite  for  tight  formation  analysis  in  320km, 
15°  orbit,  6po  =  0.0001° 
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Figure  E5. 


Figure  E6. 
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Figure  E8. 
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Figure  Ell. 


Figure  E12. 
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Figure  E14. 


Chiel 
Sail 
Sat  2 
Sat  3 
Sat  4 


phi,  (deg) 


Initial  position  in  torus  space  and  inertial  cartesian  space  of  satellite  cluster 
for  tight  formation  analysis  in  630km,  30°  orbit,  5ip  =  0.001° 


Satellite  #  1 


Satellite  #  2 


Cluster  distance  from  chief  satellite  for  tight  formation  analysis  in  630km, 
30°  orbit,  S(pQ  =  0.001° 


122 


Figure  E15. 


Cluster  drift  from  initial  separations  for  tight  formation  analysis  in  630km, 
30°  orbit,  =  0.001° 
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